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Abstract 

This  research  effort  identifies  and  models  the  nonlinear  time-variant  behavior 
exhibited  by  an  avalanche  photodiode  (APD)  array  based  laser  ranging  and  detec¬ 
tion  (LADAR)  system.  This  examination  was  prompted  by  the  failure  of  the  range 
accuracy  results,  in  test  measurements,  to  approach  the  expected  Cramer-Rao  lower 
bound  for  accuracy,  and  anomylous  signal  level  behavior  which  was  not  predicted  by 
the  system  model. 

Based  on  the  original  linear  time-invariant  (LTI)  model,  the  expected  evolution 
of  error  in  the  LADAR  signal  is  examined  sequentially  from  the  outgoing  pulse  through 
signal  digitization.  This  error  evolution  shows  the  LTI  model  does  not  contain  a 
mechanism  for  causing  the  observed  signal  deviations  or  to  explain  the  failure  of  the 
system  to  meet  the  Cramer-Rao  lower  bound.  Therefore,  a  nonlinear  time-variant 
model  is  developed  based  on  the  interactions  of  the  APDs  in  the  array  with  the 
array’s  voltage  regulator. 

In  the  refined  model,  the  sum  photo-current  for  the  entire  array  places  a  load 
on  a  voltage  regulator.  Given  a  changing  load  during  the  36  nanosecond  range  gate, 
the  voltage  regulator  cannot  maintain  its  output  at  the  reference  voltage.  The  re¬ 
sulting  reverse  bias  voltage  variations  cause  the  responsivity  of  each  APD  to  vary 
in  a  nonlinear  fashion.  Because  each  APD  in  the  array’s  responsivity  depends  upon 
the  entire  array’s  photonic  loading,  each  individual  APD’s  response  is  time  variant. 
However,  if  the  system’s  response  to  photonic  loading  is  incorporated  into  the  model 
for  each  detector,  ranging  algorithms  based  on  linear  systems  models  can  be  used  on 
the  resulting  signals. 
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Nonlinear  Time- Variant  Response 


in  an  Avalanche  Photodiode  Array  Based 
Laser  Detection  and  Ranging  System 

I.  Introduction 

1.1  Laser  Radar  Applications 

Advances  in  airpower  have  assured  the  future  of  laser  radar  on  the  modern  bat¬ 
tlefield.  In  the  absence  of  air  superiority,  an  army  will  find  that  any  vehicle,  structure, 
or  even  troop  concentration,  which  can  be  located,  can  be  destroyed.  Thus,  camou¬ 
flage  and  deception  have  taken  on  a  key  position  in  military  operations.  Keeping 
pace  with  these  developments,  modern  camouflage  extends  beyond  the  familiar  vi¬ 
sual  camouflage  patterns  which  break  up  the  outlines  of  personnel  and  equipment.  It 
has  evolved  into  full  signature  management  systems,  addressing  how  the  camouflage 
matches  its  intended  environment  in  the  visible,  infrared,  and  radar  regimes  [14]. 
Laser  Detection  and  Ranging  (LADAR)  is  a  method,  under  continuing  development, 
for  defeating  this  type  of  camouflage. 

Currently  fielded  optical  sensor  systems  primarily  rely  upon  their  spatial  res¬ 
olution  to  locate  and  identify  targets.  An  operator  examines  an  image  with  finite 
resolution  collected  by  systems  such  as  low  light  cameras,  photomultiplier  based  pas¬ 
sive  night  vision  devices,  or  forward  looking  infrared  systems  to  identify  targets.  While 
these  systems  may  provide  magnification  and  extend  the  portion  of  the  electromag¬ 
netic  spectrum  the  operator  can  view,  the  images  are  still  interpreted  by  the  operator 
for  target  recognition.  Therefore,  classical  techniques  to  deceive  operators  attempt 
to  make  the  image  of  the  object  appear  to  be  something  else  or  to  cause  insufficient 
contrast  between  the  target  and  its  background  for  detection. 

LADAR  may  be  exploited  to  overcome  both  methods  of  deception.  First,  by 
adding  range  information  to  the  image,  the  three  dimensional  shape  of  objects  can 
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be  determined  regardless  of  the  intensity  image  mapped  over  its  surface.  Even  if  the 
exact  image  of  the  terrain  behind  it  were  displayed  on  the  imaged  side  of  a  target,  its 
range  profile  will  stand  out  from  the  true  background.  If  the  material  used  to  mask  or 
break  up  the  target’s  signature  does  not  cover  it  completely,  as  is  often  the  case  with 
camouflage  netting,  there  will  be  gaps  through  which  the  target  can  be  viewed  [14]. 
LADAR  systems  capable  of  detecting  multiple  returns  can  form  multiple  layer  range 
profiles  showing  both  the  camouflage  and  what  is  underneath  it. 

To  be  effective  on  the  modern  battlefield,  a  LADAR  system  must  be  able  to 
quickly  and  accurately  make  range  measurements  on  its  entire  field  of  view  (FOV). 
As  with  all  measurement  devices,  random  and  systematic  error  is  introduced  into  the 
results.  While  typically  random  errors  can  only  be  suppressed,  systematic  errors  can 
be  removed  from  a  well-modeled  system.  The  objective  of  this  work  is  to  identify  and 
delineate  the  random  from  the  systematic  errors  in  the  range  data  generated  by  an 
experimental  FLASH  LADAR,  and  to  develop  a  model  which  can  be  used  to  remove 
the  systematic  errors. 

1.2  AFRL/SNJM  FLASH  LADAR  Specifications 

The  imaging  angle-angle-range  LADAR  operated  by  the  Electro-Optic  Com¬ 
bat  Identification  Branch  of  the  Air  Force  Research  Laboratories’  Sensors  Direc¬ 
torate  (AFRL/SNJM)  is  a  FLASH  LADAR  system  from  Advanced  Scientific  Con¬ 
cepts  (ASC)  mated  with  a  1.55-/im  Big  Sky  CFR-400  laser.  In  angle-angle-range 
format,  as  in  a  digital  camera,  a  common  optical  system  forms  an  image  on  an  array 
of  individual  detectors  or  pixels  arranged  in  a  grid  pattern.  Consequently,  each  pixel 
is  exposed  to  an  instantaneous  field  of  view  (IFOV)  spanning  some  finite  vertical  and 
horizontal  angle  within  the  FOV  of  the  entire  system. 

Range  data  is  found  by  illuminating  the  target  area  with  a  laser  pulse  and 
delaying  data  capture  to  collect  the  returning  pulse.  The  irradiance  on  each  pixel  in 
the  array  is  then  recorded  at  a  1.862-nanosecond  sampling  period  to  form  a  series  of 
nineteen  images  over  approximately  35  nanoseconds.  The  time  series  of  irradiance 
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values  recorded  by  each  pixel  is  processed  to  determine  a  range  for  the  object  or 
objects  within  its  IFOV. 

The  most  novel  component  of  the  FLASH  LADAR  system  is  an  avalanche  pho¬ 
todiode  (APD)  array  fabricated  by  Sensors  Unlimited  Inc.  which  is  mated  to  an  ASC 
Read-Out  Integrated  Circuit  (ROIC).  The  detector  array  is  128  x  128  pixels,  repre¬ 
senting  a  significant  advance  in  indium  gallium  arsenide  APD  array  fabrication  at  the 
time  of  its  development.  Due  to  the  guard  ring  structure  of  each  APD,  the  active 
detector  area  is  approximately  40- /irn  x  40-/im  with  detectors  on  100-/rm  centers.  Im¬ 
mersion  lenses  are  used  to  increase  the  effective  fill  factor  of  the  detector  array  to  near 
100  percent,  significantly  increasing  the  flux  collecting  efficiency  of  the  receiver  [17]. 

The  FLASH  LADAR  has  two  modes  for  sampling  the  period  in  which  the  laser 
pulse  actually  returns  to  the  detector.  In  stop  mode,  a  buffer  is  cyclically  filled  and 
overwritten  for  several  samples  after  any  pixel  crosses  an  irradiance  threshold.  Thus, 
the  range  to  the  target  need  not  be  guessed  to  determine  when  to  stop  overwriting 
the  buffer.  In  sular  mode,  a  software  selected  range  gate  is  selected  to  delay  write-out 
until  a  returning  pulse  from  the  selected  range  is  expected  to  arrive.  Both  modes 
result  in  an  output  record  in  which  the  first  image  written  out  is  not  necessarily  the 
first  in  the  true  time  sequence.  This  is  simply  rectified  by  reordering  the  images  after 
the  buffer  stop  time  to  be  first,  and  then  appending  any  leading  images  to  the  end  of 
the  record  [17]. 

The  sular  operating  mode  was  used  almost  exclusively  for  collecting  data  sets 
for  algorithm  development.  In  this  mode,  the  range  gate  consistently  observes  the 
same  area,  so  multiple  frames  can  be  averaged  to  reduce  noise.  The  measurement 
consistency  also  led  to  the  development  of  a  direct  method  for  removing  ambient 
light  and  fixed  pattern  noise.  This  method  consisted  of  subtracting  a  frame  of  data 
taken  in  sular  mode  with  the  laser  inactive  from  frames  taken  with  the  laser  active. 

A  system  model  was  derived  for  the  FLASH  LADAR,  incorporating  a  series  of 
assumptions,  which  are  typical  of  cameras  and  infrared  imagers.  First,  each  detector 
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was  considered  to  be  a  completely  independent  device  comprised  of  an  APD,  capacitor 
bank,  amplifier,  and  analog-to-digital  conversion  circuitry.  Second,  each  device  has 
a  characteristic  average  responsivity,  due  to  the  APD  responsivity  and  amplifier, 
and  there  is  some  variation  in  responsivity  between  devices.  Third,  each  detector 
is  capable  of  collecting  and  digitizing  the  photonic  signals  of  interest  with  negligible 
impulse  response  effects.  This  model,  for  a  single  detector,  is  shown  in  Figure  1.1. 


Target 

Scene 


Photons 


Capacitors 


♦apdHHHK]- 


10  amps/watt 


Op. 

Amp 


AID  Conversion 
&  Storage 


Figure  1.1:  Original  system  model 

1.3  Design  Impact  on  System  Performance 

The  FLASH  LADAR  system  was  designed  to  meet  a  series  of  design  require¬ 
ments.  Chief  among  these  was  a  high  a  responsivity  which  peaked  in  the  near  infrared 
region.  A  high,  for  a  developmental  device,  spatial  resolution  was  also  required. 

Imaging  arrays  of  APDs  are  significant  for  laser  radar  applications  because  they 
have  greater  than  unity  detector  gain.  The  gain  is  derived  from  one  incident  photon 
generating  many  electrons  in  the  detector  material,  and  allows  for  detector  responsiv- 
ities  as  high  as  10  Amps/ Watt  at  the  design  wavelength  of  the  APD  array  [2],  A  high 
responsivity  is  critical  for  laser  radars  examining  distant  targets,  because  very  few 
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photons  will  be  collected  by  the  imaging  optics.  Alternatively,  this  high  responsivity 
can  be  exploited  to  allow  eye-safe  experimentation  on  nearby  targets. 

While  a  high  responsivity  improves  the  signal  to  noise  ratio  of  the  LADAR 
return,  the  type  of  ranging  algorithms  which  produces  the  most  accurate  results  de¬ 
pends  strongly  on  the  system  resolution  trade  offs.  For  a  given  system,  the  highest 
possible  spatial  and  range  resolution  for  a  given  system  is  desirable,  but  any  realizable 
system  will  always  have  some  total  sample  count  limitation  per  range  image.  This 
limit  may  arise  from  repetition  rate  requirements  or  the  limitations  of  focal  plane 
array  FPA  and  ROIC  fabrication  technologies.  In  the  case  of  the  FLASH  LADAR, 
the  FPA  was  fabricated  with  128x128  detectors  and  the  ROIC  with  19  integrating 
capacitors  switching  at  a  fabrication  limit. 

Khourey  et  al.  defines  a  voxel  as  the  volume  enclosed  by  the  sample  spacing 
across  the  image  plane  and  in  range.  It  is  common  to  design  a  system  to  achieve  a 
nearly  cubic  voxel  at  a  design  range,  for  instance  one  cubic  meter  at  ten  kilometers  [7]. 
In  this  case,  with  an  f/2  lens  system,  the  sampled  voxel  is  approximately  1”  by  1” 
by  1.5’  at  a  typical  target  range  of  approximately  80  yards.  Consequently,  the  data 
collected  for  algorithm  development  has  much  more  spatial  resolution  than  range 
resolution,  due  to  sampling.  Therefore,  early  algorithmic  work  focused  on  algorithms 
capable  of  making  sub-sample  accuracy  range  determinations. 

1.4  Initial  Research  Results 

Initial  research  emphasized  target  sets  for  which  the  results  were  readily  pre¬ 
dicted  by  the  system  model.  One  of  the  first  test  was  conducted  on  a  two  surface 
target  board  set  shown  in  Figure  1.2.  Since  both  target  surfaces  are  retroreflecting, 
the  returning  signals  from  them  are  expected  to  be  delayed  and  attenuated  versions 
of  the  outgoing  pulse  shape,  with  additional  constant  valued  contributions  from  the 
ambient  light  and  dark  current.  The  return  from  any  off  target  pixel  is  expected  to 
consist  of  the  ambient  light  and  dark  current,  as  well  as  a  small  contribution  from 
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atmospheric  backscatter  of  the  laser  pulse.  An  unmodified  example  of  the  expected 
first  and  second  target  surface  returns  at  the  receiver  is  shown  in  Figure  1.3. 


Figure  1.2:  Two  surface  retroreflecting  target 

These  signals  are  conditioned  by  a  background  subtraction  scheme,  in  which  an 
ambient  light  illuminated  record  is  subtracted  from  a  pulse  illuminated  record.  An 
ambient  light  response  is  taken  by  disabling  the  laser  output  and  recording  the  system 
response  at  the  same  range  delay  as  the  pulse  illuminated  response.  Subtracting  the 
ambient  light  illuminated  signal  was  expected  to  remove  both  the  photonic  contribu¬ 
tions  due  to  ambient  light,  and  the  dark  current  response  of  the  system.  The  expected 
final  results  are  shown  in  Figure  1.4 

The  first  sub-range  sample  accuracy  algorithm  tests  on  this  target  resulted  in 
the  two  surfaces,  which  were  40”  apart,  being  measured  to  be  44”  apart.  This  result 
was  a  significantly  less  accurate  measurement  than  was  expected  for  the  algorithm 
being  tested,  which  was  based  on  a  linear  time-invariant  system  model  and  postulated 
return  pulse  shape  [1] .  While  examining  these  results,  it  was  noted  that  the  recorded 
signals  did  not  match  the  signals  predicted  by  the  system  model.  The  recorded  returns 
are  displayed  in  Figure  1.5. 
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Figure  1.3:  a.  First  surface  return  b.  Second  surface  return  c.  Off  target 

While  it  appears  that  there  are  two  returns  in  each  of  pulse  illuminated  returns, 
this  cannot  be  the  case.  The  targets  in  the  IFOV  of  the  pixels  shown  are  retroreflect- 
ing.  This  indicates  that  the  observed  behavior  may  be  some  form  of  system  response 
ringing.  Upon  close  inspection,  it  can  also  be  seen  that  the  appearant  pulse  widths 
are  different,  between  the  first  and  second  surfaces.  Once  again,  this  cannot  be  the 
case  for  a  pair  of  retroreflecting  surfaces.  Finally,  the  behavior  of  the  off  target  pixel 
deviates  severely  from  the  expected  system  response.  Any  response  due  to  atmo¬ 
spheric  backscatter  must  begin  as  the  laser  pulse  leaves  the  aperture,  and  will  be  a 
nearly  linear  response  when  measured  at  a  significant  distant  from  the  aperture,  over 
the  range  gate. 

An  additional  data  set  under  consideration  consists  of  the  same  target  board 
set,  in  a  much  more  range  diverse  scene.  In  this  case,  the  first  surface  of  the  target 
represents  the  start  of  a  broader  set  of  range  returns  including  the  second  surface, 
personnel,  and  foliage.  Figure  1.6  subplot  a.  shows  the  pulse  shape  recorded  for  the 
first  surface  of  the  target  board.  This  pulse  exhibits  a  broader  shape  than  expected, 
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Figure  1.4:  a.  First  surface  return 


b 


.  Second  surface  return  c.  Off  target 


and  what  may  be  the  start  of  the  anomylous  ringing  behavior.  In  the  off  target  pixel 
response  in  Figure  1.6  subplot  b.,  a  high  negative  correlation  can  be  seen  between  the 
value  after  background  subtraction  and  the  pulse  shape  from  the  first  surface.  This 
response  is  distinctly  different  from  the  ambient  response  value,  after  dark  current 
subtraction,  shown  in  Figure  1.6  subplot  c. 

This  anomylous  behavior  has  several  significant  implications  on  the  military 
application  of  this  LADAR  system.  The  observed  behavior  makes  multiple  object 
ranging  inside  an  IFOV  impossible.  Any  multiple  surface  ranging  algorithm  would 
assign  two  distinct  ranges  to  the  recorded  data  in  Figure  1.5,  even  though  a  second 
surface  is  known  to  be  absent.  In  general,  ranging  algorithms  would  also  assign  a  range 
to  the  off  target  pixels,  most  likely  near  the  twelfth  sample  in  Figure  1.5.  However, 
the  effects  of  these  behaviors  on  the  accuracy  of  single  surface  ranging  was  uncertain. 

In  an  effort  to  develop  a  better  expectation  of  the  system  and  algorithm’s  an¬ 
alytical  measurement  accuracy,  the  Cramer-Rao  lower  bound  was  calculated  for  the 
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Figure  1.5:  Pulsed  values  minus  background  a.  First  surface  return  b.  Second 

surface  return  c.  Off  target  return 

sub-range  sample  accuracy  algorithm  under  consideration,  given  a  single-surface  test 
case.  For  this  case,  the  bound  was  found  to  be  approximately  3/5”.  However,  this 
result  is  for  an  unbiased  estimator  and  the  ranging  algorithm  does  not  meet  this  re¬ 
quirement  for  returns  exhibiting  the  observed  anomylous  behavior  [16].  In  practice, 
the  FLASH  LADAR  range  accuracy  for  the  algorithm  and  test  case  was  found  to 
be  2.5”.  Therefore,  the  next  major  research  objective  is  to  identify  the  source  of  and 
model  the  anomylous  signal  behavior,  and  implement  the  resulting  model  into  ranging 
algorithms. 


1.5  Research  Scope 

The  primary  objective  of  this  work  is  to  investigate  and  model  the  anomylous 
signal  behavior  of  the  FLASH  LADAR  under  optical  loading.  A  refined  system  level 
model  is  proposed,  based  upon  key  system  elements  for  which  partial  models  or  de¬ 
scriptions  are  available.  Due  to  the  lack  of  an  appropriately  controlled  test  case,  only 
a  forward  model  is  developed.  Although  the  detailed  analysis  to  follow  is  specific  to 
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Figure  1.6:  a.  First  surface  after  background  subtraction  b.  Off  target  pixel  after 
background  subtraction  c.  off  target  ambient  response  minus  dark  response 

the  FLASH  LADAR,  the  general  principles  of  the  model  are  applicable  to  any  APD 
array  based  LADAR,  of  similar  design. 

1.6  Research  Methodology 

The  type  and  scope  of  investigations  which  can  be  carried  out  to  address  this 
problem  are,  by  necessity,  limited.  The  functionality  of  the  LADAR  system  cannot 
be  compromised,  so  partial  disassembly  for  component  measurements  is  not  possible. 
There  is  also  a  limited  amount  of  published  information  regarding  detailed  system 
specifications,  due  to  intellectual  property  and  security  requirements.  Consequently, 
the  system  is  treated  in  a  gray-box  manner,  with  limited  known  parameters. 

The  analytical  approach  which  follows  was  used  to  identify  the  deviations  be¬ 
tween  the  modeled  and  recorded  system  responses,  and  develop  a  refined  model.  First, 
the  original  system’s  error  model  is  presented  following  the  laser  pulse  from  trans¬ 
mission  to  and  from  the  target,  through  the  optical  system,  into  the  detection  and 
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digitization  process,  and  finally  through  range  interrogation  by  various  algorithms. 
This  signal  error  evolution  accounts  for  deviations  due  to  known  sources,  in  their 
generally  accepted  forms.  Then,  a  hypothetical  system  model  is  proposed  to  explain 
the  anomylous  behavior.  Beginning  with  the  digitized  output  of  the  system,  the  de¬ 
tection  process  is  reversed  to  examine  the  system  behavior,  and  develop  a  predictive 
forward  model  for  the  system.  The  forward  model  is  completed  by  solving  a  system  of 
linear  constant-coefficient  difference  equations  to  identify  the  remaining  linear  system 
behavior.  Finally,  the  proposed  forward  model  is  compared  to  a  recorded  case  for 
validation. 
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II.  Background 


2. 1  Introduction 

The  objective  of  the  FLASH  LADAR  system  is  to  assign  a  range  to  each  pixel  by 
evaluating  a  recorded  pulse  return  with  a  ranging  algorithm.  It  was  presumed,  in  the 
original  model,  that  the  LADAR  system’s  response  was  both  linear  time-invariant  and 
had  negligible  effects  on  the  recorded  pulse  shape.  The  performance  of  any  ranging 
algorithm  is  degraded  by  random  noise,  and  all  ranging  algorithms  incorporate  a 
series  of  assumptions  about  the  signal  being  processed.  If  the  noise  level  is  too  high 
or  the  assumptions  are  in  error,  the  Cramer- Rao  limit  for  ranging  accuracy  cannot  be 
achieved.  Therefore,  the  expected  evolution  of  error  in  the  FLASH  LADAR  model 
is  examined  to  account  for  known  systematic  and  random  effects.  Then,  a  series 
of  common  ranging  algorithms  and  their  associated  assumptions  are  presented  and 
compared  to  the  behavior  of  the  recorded  signals. 

2.2  Signal  Error  Evolution 

The  ideal  return  signal  for  the  FLASH  LADAR  is  the  convolution  of  a  linearly 
scaled  version  of  a  known  outgoing  laser  pulse  shape  with  a  single  retro-reflecting 
surface  in  each  pixel’s  IFOV  [18].  The  resulting  return  is  an  attenuated  copy  of  the 
outgoing  laser  pulse  recorded  some  time  after  the  pulse  was  transmitted.  However, 
numerous  factors  cause  the  recorded  return  to  deviate  from  this  ideal  state:  random 
error  is  introduced  by  the  lasing  process;  the  atmosphere  between  the  FLASH  LADAR 
and  target  surfaces;  the  photonic  detection  process,  and  the  digitization  process. 
Systemic  errors  are  also  introduced  by  the  range  profile  of  the  target  within  each  pixel’s 
IFOV,  the  system’s  optical  train,  the  detector  array,  and  the  digitization  process. 
These  sources  of  error  will  be  examined  following  the  system’s  signal  path  from  the 
outgoing  laser  pulse  to  the  reported  digital  output. 

2.2.1  Transmitted  Pulse.  A  well-characterized  laser  source  is  a  critical  first 
component  of  any  LADAR  system,  because  all  ranging  methods  require  some  knowl- 
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edge  of  the  transmitted  pulse  shape.  While  well-behaved  pulses  can  be  consistently 
generated,  there  will  always  be  some  random  error  between  the  designed  and  the 
actual  individual  outgoing  pulse  shape.  When  ranging  algorithms  are  applied,  any 
deviation  from  the  expected  pulse  shape  in  a  LADAR  return  will  eventually  mani¬ 
fest  itself  as  a  range  estimation  error.  Even  the  mathematically  ideal  returning  pulse 
shape  can  introduce  error  in  deconvolution,  if  more  than  one  range  is  represented  in 
a  pixel’s  IFOV  [18]. 

The  transmitted  pulse  must  be  of  sufficient  power  that  a  detectable  number  of 
photons  will  be  scattered  back  towards  the  LADAR  from  the  target  being  illuminated. 
Additional  power,  beyond  the  device  minimum  for  detetion,  is  required  because  the 
signal  entering  the  detector  will  not  consist  solely  of  the  laser  pulse  energy  reflected 
from  the  target,  ft  will  also  include  in-band  photons  from  other  sources  such  as  self 
emission  or  solar  reflection  [9].  The  pulse  must  contain  sufficient  power  to  achieve  at 
least  a  signal-to-noise  ratio  of  one  with  these  sources,  upon  arrival  at  the  detector. 

The  atmosphere,  through  which  the  laser  pulse  travels,  affects  the  laser  pulse 
traveling  both  to  and  from  the  target.  Scattering  and  absorption  in  the  atmosphere 
attenuate  the  returning  pulse  by  reducing  the  number  of  photons  which  return  to  the 
FLASH  LADAR.  Speckle  noise,  which  arises  from  the  coherence  of  the  laser  pulse,  also 
causes  error  in  the  returning  signal.  The  end  result  of  speckle  noise  is  to  cause  random 
deviations  in  the  flux  distribution  collected  at  the  receiver,  due  to  interference  [4], 

2.2.2  Optical  Train.  The  first  systematic  error  introduced  by  the  receiver  is 
due  to  the  optical  system’s  point  spread  function  (PSF).  For  perfect  optics,  the  PSF 
is  governed  by  diffraction  through  the  finite  aperture  of  the  system’s  entrance  pupil. 
The  image  formed  by  the  optical  system  is  the  convolution  of  the  PSF  with  the  flux 
distribution  on  the  system’s  entrance  pupil.  The  normal  result,  for  physical  systems, 
is  to  introduce  light  from  a  point  source  onto  a  finite  area. 

Assuming  diffraction  limited  performance,  the  f/2  lens  system  used  is  expected 
to  yield  a  diffraction  blur  spot  diameter  ddifr  of  7.564-/im  according  to  the  relationship 


13 


expressed  in  equation  (2.1)  [3].  For  the  discretely  sampled  image  plane  of  the  LADAR, 
the  key  deviation  arises  from  a  portion  of  the  irradiance  which  was  expected  to  fall 
on  a  certain  pixel  instead  falling  on  multiple  nearby  pixels.  Therefore,  if  uncorrected, 
the  irradiance  record  on  each  pixel  will  contain  information  from  outside  its  1FOV  [5]. 

ddiff  =  2.44  *  A  *  f/#  (2.1) 

Any  realizable  lens  system  also  introduces  aberrations  in  the  image.  The  laser 
source  means  that  the  light  of  interest  entering  the  system  is  quasi-monochromatic, 
and  the  primary  aberrations  will  be  the  Seidel  aberrations  [6].  For  the  purposes  of 
ranging,  the  most  harmful  of  effect  spherical  aberration,  coma  and  astigmatism  is 
to  blur  the  image  by  introducing  light  from  outside  each  pixel’s  1FOV.  Petzval  field 
curvature  and  distortion  also  redistribute  the  light  from  the  target.  If  the  image  is 
properly  sampled,  these  aberrations  spatially  misplace  the  resulting  range  from  its  true 
location  on  the  object,  but  may  not  introduce  ambiguity  in  the  ranging  process  [6]. 

The  average  PSF  can  be  accurately  measured,  and  the  aberrations  minimized 
and  measured  for  an  optical  train.  This  is  particularly  straight  forward  for  sys¬ 
tems  such  as  LADARs,  which  are  nearly  monochromatic  [6].  For  well-characterized 
systems,  these  effects  can  be  deconvolved  from  the  data  if  the  minimum  sampling 
requirements  are  met  and  random  noise  can  be  addressed  effectively  [5]. 

2.2.3  Photonic  and  Detector  Noise.  Photonic  detection  is  a  random  arrival 
process,  which  is  well  studied  for  the  system  under  consideration  [4].  The  pulse 
returns  from  a  coherently  illuminated  scene  onto  a  detector  with  a  sampling  area 
which  is  much  larger  than  a  coherence  area  on  the  detector.  Therefore,  the  arrival 
process  takes  on  a  negative  binomial  form  dependent  upon  the  average  photon  count 
and  speckle  parameter  [16].  The  system  operates  at  f/2,  with  a  center  wavelength  of 
1.55-/im,  and  a  detector  size  of  100-pm  x  100-pm.  This  yields  a  speckle  cell  width  of 
approximately  3- /mi,  and  a  detector  area  to  speckle  cell  area  ratio  of  178  to  1.  This 
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ratio  is  the  speckle  parameter.  For  speckle  parameters  much  greater  than  one,  such 
as  this  case,  the  arrival  distribution  approaches  a  Poisson  distribution  [16]. 

The  detectors  in  the  FLASH  LADAR  are  avalanche  photodiodes.  APDs  have 
a  non-unity  internal  gain  term  which  multiplies  shot  and  photonic  noise  along  with 
the  signal,  and  introduces  its  own  excess  noise  term  as  defined  in  equation  (2.2). 
This  excess  noise  term,  F.  stems  from  the  statistical  nature  of  the  gain  term,  and  is 
mathematically  expressed  as  a  function  of  the  average  gain,  G,  the  electron  impact 
ionization  coefficient,  ot,  and  the  hole  impact  ionization  coefficient,  (3  [3].  The  excess 
noise  factor  typically  begins  to  dominate  other  noise  sources,  including  read-out  noise, 
above  gains  of  twenty.  It  is  not  expected  to  be  the  dominant  noise  source  for  this 
system,  which  has  a  gain  of  approximately  ten. 


(2.2) 


Ideally,  if  each  detector  in  the  array  were  exposed  to  the  same  incident  flux, 
it  would  report  the  same  averaged  photo-electron  count.  In  reality,  each  APD  in 
the  array  has  its  own  average  gain  term  and  dark  current.  For  constant  operating 
conditions,  these  terms  are  dependent  upon  the  physical  parameters  of  each  device. 

When  exposed  to  a  common  input  flux  and  allowed  to  reach  a  steady-state 
value,  the  multiplicative  effect  of  each  device’s  gain  and  the  additive  effect  of  its 
dark  current  can  be  divided  and  subtracted  out  to  produce  a  common  photo-electron 
count.  The  terms  necessary  for  this  normalization  can  then  be  recorded  for  each  pixel, 
as  a  correction  factor.  This  type  of  correction  is  typically  termed  non-uniformity 
correction  (NUC)  and  is  conducted  on  the  final  digital  output  of  the  device  to  remove 
the  combined  deviations  of  the  detector,  amplifier,  and  sampling  circuitry  [9]. 

2.2.4  Sampling  in  Space  and  Range.  The  detector  of  an  angle-angle-range 
LADAR  must  sample  intensities  in  both  dimensions  of  the  image  plane  and  in  range. 
The  Nyquist  sampling  requirement  must  be  met  in  each  dimension  to  form  an  un- 
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aliased  signal  at  the  designed  resolution.  The  factors  determining  the  Nyquist  sample 
rate  are  quite  different  between  the  image  plane  and  in  range. 

Across  the  image  plane,  the  signal  is  truly  band-limited  due  to  the  system’s 
entrance  pupil,  and  the  signal  can  be  sampled  at  a  finite  rate  with  no  aliasing.  For 
a  single  wavelength,  the  highest  spatial  frequency,  fc,  which  is  passed  by  the  optical 
train  can  be  calculated  via  equation  (2.3).  In  this  equation,  u  is  the  effective  aperture 
radius,  z,  is  the  distance  from  the  aperture  to  the  focal  plane,  and  A  is  the  wavelength 
under  consideration. 

The  FLASH  LADAR  uses  an  f/2  lens  with  a  10-cm  entrance  pupil.  Therefore, 
the  highest  spatial  frequency  which  is  passed  by  the  optical  system  is  fc  =  1.613  *  105 
mT1.  This  requires  a  Nyquist  sampling  rate  of  2*/c  which  corresponds  to  a  spacing  of 
3.1-/im.  The  system  is  sampled  at  100- /rm  spacing,  and  is  thus  spatially  undersampled. 
Consequently,  all  of  the  information  recorded  by  the  detector  array  exhibits  spatial 
aliasing. 


(2-3) 

In  range,  the  pulse  shape  defines  the  Nyquist  sampling  requirement  and  thus  the 
signal  is  never  truly  band-limited.  As  always  when  sampling  a  non-band  limited  sig¬ 
nal,  some  noise  is  introduced  due  to  aliasing.  Typically,  the  aliasing  noise  is  mitigated 
by  low  pass  filtering  the  signal  prior  to  digitization.  In  this  case,  the  smoothly  vary¬ 
ing  pulse  shape  provides  significant  high-frequency  attenuation,  reducing  the  aliasing 
effect.  Additionally,  any  realizable  optical  detector  applies  a  low-pass  filtering  oper¬ 
ation  to  the  incident  signal,  due  to  the  finite  detector  integration  time  [5].  For  this 
system,  the  actual  integration  time  of  the  detector  is  not  known  explicitly,  although 
it  is  necessarily  shorter  than  the  switching  time  of  1.862  nanoseconds.  Based  on  these 
factors,  the  noise  due  to  aliasing  in  range  is  assumed  to  be  negligible  when  compared 
to  other  noise  sources  for  this  system. 
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2.2.5  Read-Out  and  Digitization.  In  the  temporal  sampling  process,  each 
detector  sequentially  charges  twenty  integrating  capacitors.  At  the  completion  of  the 
overwriting  cycle,  the  charge  in  the  capacitors  is  read  out  via  an  amplifier  and  under¬ 
goes  an  analog-to-digital  conversion  process.  This  process  introduces  read-out  noise 
based  on  charge  deviation,  individual  amplifier  gain  variation,  individual  amplifier 
additive  offsets,  injected  amplifier  noise,  and  the  rounding  process  associated  with 
digitization.  The  random  error  introduced  during  read  out  is  typically  modeled  as 
Gaussian  White  Additive  Noise,  GWAN,  while  the  variation  between  sets  of  detec¬ 
tor  circuitry  is  represented  by  an  array  of  multiplicative  and  additive  terms.  Both 
of  these  types  of  error  enter  the  final  digital  count  stored  by  the  LADAR,  which  is 
significantly  rescaled  in  comparison  to  a  true  photo-electron  count  unless  the  system 
is  calibrated  against  a  known  photonic  source. 

2.3  Range  Determination 

An  ever-increasing  number  of  methods  for  deriving  range  estimates  from  LADAR 
returns  are  available.  Most  of  these  methods  are  fundamentally  deconvolution  meth¬ 
ods  with  different  means  of  addressing  uncertainty  in  the  target  shape  and  noise.  Four 
common  methods  which  will  be  examined  include  peak  detection,  matched  filtering, 
Richardson- Lucy  deconvolution,  and  an  estimation-theory-based  method  [18]. 

2.3.1  Peak  Detection.  The  simplest  method  of  range  estimation  is  peak 
detection.  Peak  detection  algorithms  are  best  matched  with  pulses  shorter  than  the 
sampling  time  of  the  system,  but  can  be  applied  to  longer  pulsed  systems.  This 
method’s  range  accuracy  is  directly  defined  by  the  clock  sampling  time,  since  the  peak 
irradiance  value  will  always  be  localized  in  a  single  time  increment.  If  the  returning 
pulse  is  temporarily  recorded,  a  clock  is  started  when  the  pulse  is  fired.  Then,  when 
the  irradiance  on  the  detector  crosses  a  preset  threshold,  the  clock  is  stopped.  The 
resulting  round  trip  time  is  converted  to  a  round  trip  distance,  assuming  the  speed  of 
light  in  a  vacuum.  If  data  is  recorded  shortly  before  and  after  the  stop  time,  multiple 
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threshold  crossings  can  be  evaluated.  In  this  way,  multiple  reflections  in  the  field  of 
view  can  be  ranged  rather  than  only  the  first  [18]. 

Peak  detection  provides  consistent  results  in  the  presence  of  noise,  but  provides 
the  least  accurate  range  estimation  of  all  the  methods  compared  in  this  section.  Using 
peak  detection,  the  FLASH  LADAR  system  has  a  range  uncertainty  of  approximately 
1.5’,  which  is  acceptable  for  weapons  ranging,  but  does  not  provide  the  detail  necessary 
for  LADAR  based  target  identification  [16]. 

2.3.2  Matched  Filtering.  In  matched  filter  detection,  which  is  commonly 
used  in  radar  processing,  a  noisy  or  deformed  input  signal  is  processed  by  a  filter 
whose  impulse  response  is  the  time  reverse  of  the  signal  to  be  detected.  The  peak 
of  the  resulting  filtered  waveform  represents  the  point  in  time  at  which  the  highest 
correlation  was  achieved  [10].  The  mathematical  implementation  of  this  system  is 
very  fast,  but  multiple  returns  in  the  same  signal  and  other  deviations  from  the  ex¬ 
pected  return  shape  significantly  impact  the  accuracy  of  this  method.  Any  additional 
signal  content  shifts  the  position  of  highest  correlation,  and  thus  the  range  estimate, 
towards  the  deviation  signal’s  range.  Random  noise  is  expressed  by  correlation  peak 
broadening  and  results  in  ambiguity  in  the  range  estimate. 

2.3.3  Richards  on- Lucy  Deconvolution.  Many  current  deconvolution  algo¬ 
rithms  are  based  on  the  iterative  Richardson-Lucy  algorithm.  As  such,  deconvolution 
begins  by  numerically  estimating  the  probability  density  function  of  an  observed  wave¬ 
form,  o(q,r,k),  using  each  sample  from  each  pixel  position  defined  by  q  and  r  with 
a  range  return  centered  at  k.  The  intensity  value  in  each  time  increment  d(q,r,k) 
is  assumed  to  be  a  realization  of  a  random  variable,  D(q,r,k),  with  a  probability 
density  function  based  on  a  noise  model  appropriate  for  the  system.  Bayes  Rule  indi¬ 
cates  the  relationship  between  the  as-yet- unknown  probability  of  the  actual  waveform 
conditioned  on  the  observed  waveform,  the  probability  of  the  observed  waveform,  and 
the  probability  of  the  observed  waveform  conditioned  on  the  actual  waveform.  This 
relationship  is  shown  in  Equation  (2.4). 
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A  log  likelihood  estimator  for  o(q,r,  k )  then  takes  the  form  of  Equation  (2.5). 


L(o)  =  ln(fD\0(d\  o))  (2.5) 

The  likelihood  estimator  is  then  expanded  into  terms  determined  by  the  system 
model,  differentiated,  and  set  equal  to  zero.  An  iterative  updating  solution  is  used 
to  solve  the  resulting  intractable  equation  for  o(q,r,knew )  [15].  The  solution  of  the 
equation,  o(q,r,  knew),  is  a  representation  of  a  noisy  range  return  for  the  pixel  of 
interest  at  a  given  range,  knew.  If  the  estimated  noise  probability  distribution  function 
does  not  represent  the  actual  noise  distribution  properly,  this  method  may  converge 
on  a  false  range  or  fail  to  converge  at  all.  This  method  responds  to  bias  deviations  in 
a  similar  manner  to  matched  filter  detection. 

2.3.4  Range  Gain  and  Bias  Estimation.  Estimation  theory  methods  at¬ 
tempt  to  directly  fit  a  scaled  and  shifted  version  of  the  output  pulse  to  the  recorded 
data  and  iteratively  locate  the  pulse  position  giving  the  least-mean-squared  error 
(MSE).  This  type  of  algorithm  requires  searching  through  all  pulse  locations  under 
consideration  and  is  more  time  intensive  than  other  implementations,  but  can  esti¬ 
mate  object  ranges  outside  the  sampled  window  and  multiple  return  locations.  The 
range-gain-bias  (RGB)  estimation  method  derived  by  Burris  [1],  uses  a  three  term 
model  based  on  an  inverted  parabola  pulse  shape  with  an  initially  unknown  gain, 
location,  or  bias  level.  The  pulse  model,  I(tk,xn,ym),  is  shown  in  Equation  (2.6). 

The  gain  term,  G(xn,ym),  scales  the  pulse-shape  model,  which  is  truncated  by 
a  shifted  rectangle  function,  rect(f(tk)),  with  a  pulse  half-width,  pw,  in  seconds  mul¬ 
tiplied  by  the  speed  of  light,  c,  to  preclude  negative  pulse  values.  Noise  in  each  voxel 
is  modeled  by  the  q(tk,xn,ym )  function.  An  average  bias,  B(xn,ym),  is  estimated  for 
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each  test  location  by  using  the  mean  of  all  of  the  sampled  data  outside  the  postulated 
pulse.  Both  the  gain  and  bias  terms  are  solved  using  an  estimated  range 

I{tk,  xn,ym)  =  G(xn,ym)  (l  -  reci  ... 

•••  T  Urn )  “1“  Q{tki  Urn)  (2-6) 

Bias  level  averaging  makes  this  method  tolerant  of  signal  deviations  outside  the 
returning  pulse,  so  long  as  the  actual  return  pulse  has  higher  amplitude  than  any 
deviation  [1],  If  the  signal  deviates  during  the  pulse,  the  shape  of  the  pulse  and  the 
resulting  range  estimate  will  be  skewed.  Random  noise  introduces  ambiguity  in  the 
pulse  shape,  and  thus  the  range  estimate. 

2.4  Signal  Error  and  Ranging  Summary 

Throughout  the  course  of  its  transit,  the  LADAR  signal  has  deviated  signifi¬ 
cantly  from  its  desired  form.  Random  error  has  entered  from  the  laser  pulse  generation 
process,  atmospheric  speckle,  photonic  arrival,  APD  gain  noise,  and  read-out  noise. 
The  signal  has  also  been  systematically  deformed  by  the  target’s  range  profile  in  the 
IFOV,  the  optical  system’s  PSF  and  aberrations,  detector  non-uniformity,  and  read¬ 
out  amplifier  non-uniformity.  These  error  sources  are  intrinsic  to  the  device,  and  must 
be  dealt  with  by  the  ranging  algorithms  to  produce  an  accurate  range  estimate. 

Ranging  methods  including  peak  detection,  matched  filtering,  Richardson-Lucy 
deconvolution,  and  a  Range-Gain-Bias  estimation  method  were  examined.  Since  all 
four  methods  are  fundamentally  deconvolution  methods,  a  linear  time  invariant  (LTI) 
model  for  the  recorded  signals  was  intrinsically  assumed. 

In  their  generally  accepted  forms,  none  of  the  error  sources  or  ranging  models 
explain  the  anomylous  signal  behavior  observed.  However,  the  signal  flow  and  LTI 
requirements  for  accurate  ranging  provide  an  insight  into  a  refined  system  model 
which  can  explain  the  anomylous  behavior  and  ranging  errors. 
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III.  System  Model 


3.1  Model  Hypothesis 

Based  on  observations  from  Figure  1.5  and  Figure  1.6  and  the  results  of  the  evo¬ 
lution  of  error  analysis,  all  of  the  assumptions  indicated  by  the  original  system  model 
may  be  erronious.  First,  the  ringing  response  and  pulse  width  variation  observed  in 
Figure  1.5  indicate  a  non-negligible  detector  impulse  response.  For  the  original  model, 
a  non-negligible  impulse  response  indicates  that  the  APD  and  amplifier  responsivi- 
ties  are  not  constant  valued  between  samples  in  the  record.  Thus,  the  required  NUC 
coefficients  are  effectively  changing  during  the  record  and  introducing  a  systematic 
error.  Finally,  the  correlation  between  devices  which  have  no  common  photonic  input 
means  that  the  detectors  in  the  array  are  not  independent.  Based  on  this  hypothesis, 
a  refined  model  of  the  system  is  presented. 

3.2  Proposed  System  Model 

The  following  system  model  is  presented  to  indicate  the  relationship  between  the 
key  components  of  the  system,  as  they  relate  to  the  APD  gain.  The  model  provides 
a  framework  for  determining  the  system  behavior,  which  results  in  gain  variation.  It 
also  assigns  the  time  dependence  and  nonlinearity  exhibited  by  the  system  to  distinct 
components,  so  these  factors  can  be  removed  and  the  residual  system  identification 
problem  solved. 

A  proposed  functional  model  of  the  system  can  be  seen  in  Figure  3.1.  The 
photons  incident  on  each  pixel  are  multiplied  by  a  gain  term  to  yield  a  number  of 
electrons,  which  flow  as  current  to  the  integrating  capacitors  of  the  ROIC  [17].  The 
sum  of  each  pixel’s  current  demand  is  fed  back  to  the  voltage  regulator  to  produce  the 
bias  voltage  which  drives  the  APD  gain  equation  shown  in  Equation  (3.1).  The  charge 
stored  in  the  integrating  capacitors  is  resistively  converted  to  voltage  and  passed  to 
an  analog-to-digital  converter  to  produce  the  digital  photo-counts  which  are  stored  in 
the  system’s  memory. 
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The  configuration  of  these  components  suggests  that  a  parameterized  model  of 
the  system  could  be  developed,  but  some  information  is  missing.  While  the  voltage 
regulator  is  expected  to  be  a  closed  loop  control  device  of  the  proportional  integral 
derivative  (PID)  type,  its  system  parameters  are  unknown.  The  network  through 
which  the  current  demand  is  converted  back  to  a  reference  voltage  for  the  voltage 
regulator  is  also  unknown,  but  assumed  to  be  a  resistive-capacitive  (RC)  rather  than 
purely  resistive  network  [8].  Additional  uncertainty  is  added  by  the  RC  network  of  the 
integrating  capacitors  and  analog-to-digital  converters,  since  the  only  data  available 
is  the  output  recorded  by  the  system  after  analog-to-digital  conversion.  Models  for 
the  APD  gain,  voltage  regulation  process,  and  ROIC  are  presented  in  greater  detail 
in  sections  3.3,  3.4,  and  3.5. 


Figure  3.1:  Proposed  system  model  diagram 
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3.3  APD  Model  Details 


Fundamentally,  an  InGaAs  APD  is  a  highly  reverse  biased  PN  junction,  which 
is  approximately  thirty  percent  quantum  efficient.  When  a  photon,  of  the  correct 
wavelength,  is  incident  on  the  junction  an  electron-hole  pair  is  excited  out  of  the 
depletion  region.  The  strong  reverse  bias  causes  these  carriers  to  accelerate  and  gain 
enough  kinetic  energy  to  excite  new  carriers  by  impact  ionization  as  they  leave  the 
depletion  region.  These  additional  carriers,  which  are  generated  from  one  photon, 
constitute  the  APD  gain.  The  number  of  carriers  generated  depends  on  the  physical 
parameters  of  the  device  affecting  the  carrier  extinction  ratios  and  on  the  reverse  bias 
voltage  applied  across  the  diode  [3].  The  empirically  derived  Miller’s  Equation  is  used 
to  directly  relate  the  APD  gain  to  the  applied  reverse  bias  voltage,  . 

Miller’s  Equation,  shown  in  equation  (3.1),  indicates  that  the  average  gain,  G,  is 
a  function  of  the  reverse  breakdown  voltage,  Vbd,  the  reverse  biasing  voltage,  Veias, 
and  an  empirical  scaling  factor,  p.  The  reverse  breakdown  voltage  is  a  physical 
parameter  of  the  device,  although  it  can  vary  with  temperature  [19].  The  scaling 
factor  is  found  by  solving  equation  (3.1)  for  p  using  known  values  for  the  voltages 
and  gain.  In  this  case,  equation  (3.1)  yields  a  scaling  factor  value  of  5.48.  The  results 
of  equation  (3.2)  can  then  be  used  to  determine  the  gain  resulting  from  a  given  reverse 
bias  and  breakdown  voltage. 
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3-4  Voltage  Regulator  Model  Details 

The  reverse  bias  voltage  is  provided  to  all  of  the  APDs  in  the  array  by  a  single 
fast  transient  response  voltage  regulator.  This  device  can  be  simply  viewed  as  a 


current  controlled  voltage  device.  It  provides  a  stable  output  voltage,  when  a  stable 


current  demand  is  placed  on  it.  However,  when  the  current  demand  varies,  the  device 
has  a  transient  response  during  which  the  voltage  will  not  be  maintained  at  its  set 
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level  [12].  For  the  FLASH  LADAR,  the  current  demand  on  the  voltage  regulator 
varies  with  the  photocurrent  of  the  entire  array,  since  all  of  the  APDs  are  biased  by 
the  same  regulator. 

The  following  example  is  given  as  justification  for  modeling  a  variable  voltage 
regulator  output.  The  published  dark  current  for  one  pixel  of  a  developmental  version 
of  the  device  is  14.4-nA  [2].  Under  constant  dark  conditions,  the  current  load  required 
by  the  entire  array  would  be  16384  times  this  amount,  or  0.236  mA.  If  a  scene  were 
illuminated  by  a  pulse  from  the  specified  Big  Sky  CFR-400  laser,  which  outputs  7  mJ 
of  energy  per  pulse,  the  current  demand  must  increase  during  the  pulse  return. 

The  radiometric  calculations  required  to  determine  the  photonic  input  to  the 
system  for  a  test  board  80-rn  away  progress  in  the  following  manner.  The  vertical 
and  horizontal  half  angle  FOV  for  an  f/2  lens,  with  a  6.4-mm  wide  detector,  is  32- 
milliradians.  At  80-m  this  corresponds  to  viewing  an  area  5.122-m  x  5.122-m.  In  the 
data  set  available  for  analysis,  the  target  board  is  centered  in  the  FOV  and  covers  a 
2.43-rn  x  2.43-rn  square  area. 

Presuming  the  laser  was  focused  such  that  the  designed  FOV  was  covered  by  the 
first  two  standard  deviations  of  the  pulse’s  spatial  Gaussian  distribution,  43-percent 
of  the  transmitted  flux  is  incident  on  the  target  board.  The  white  board  is  presumed 
nominally  10-percent  reflective  and  Lambertian,  so  the  10-crn  collecting  optic  of  the 
receiver  collects  only  6.71  *  10“ '-percent  of  the  transmitted  pulse  or  4.7-nJ.  For  this 
example,  the  APD  gain  is  temporarily  assumed  to  be  constant  at  a  value  of  10.  Using 
a  5-range  sample  full  width  half-max  for  a  7-range  sample  pulse  indicates  that  the 
central  5-range  samples  of  the  pulse  contains  52  percent  of  the  temporal  pulse  energy. 
Therefore,  the  mean  current  demand  on  the  entire  array,  during  the  5-range  sample 
central  period  is  984-mA,  assuming  30-percent  quantum  efficiency.  Therefore,  during 
this  pulse  return,  the  current  demand  increased  by  a  factor  of  4169  over  the  dark 
response  level. 
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The  specific  voltage  regulator  could  not  be  identified  or  removed  for  testing,  but 
it  is  unlikely  that  device  will  be  able  to  compensate  for  a  984-rnA  demand  swing  in 
5-ns,  without  deviation.  Therefore,  the  supplied  reverse  bias  voltage  will  vary  under 
optical  loading.  The  time  scale  of  the  sampling  operation  ensures  that  the  transient 
responses  will  be  present  in  the  data,  even  though  the  voltage  regulator  attempts  to 
maintain  a  constant  reverse  bias  voltage. 

3.5  ROIC  Model  Details 

Each  APD  in  the  array  is  multiplexed  to  a  set  of  twenty  integrating  capacitors. 
During  data  collection,  photo-current  and  noise  current  from  the  APD  sequentially 
charges  these  capacitors  for  a  period  less  than  the  1.862-ns  switching  time.  Upon 
demand,  the  capacitor  charges  are  serially  recorded  through  an  amplifier  and  analog- 
to-digital  converter  to  system  memory  as  a  digital  count.  The  series  of  digital  counts 
recorded  from  the  integrating  capacitors  represents  the  irradiance  on  that  APD  at 
the  time  of  recording.  The  spatial  composite  of  all  the  pixels’  records  form  a  frame 
of  data. 

For  each  frame  of  sular  mode  data,  a  laser  pulse  is  triggered  and  a  detector 
at  the  laser’s  aperture  monitors  the  actual  pulse  output  from  the  laser.  When  the 
output  is  detected,  the  LADAR’s  internal  clock  begins  cycling  and  the  detector  array 
is  activated  and  begins  cyclically  over-writing  its  capacitor  charges.  Once  a  time  delay 
corresponding  to  the  selected  range  has  elapsed,  the  detector  array  is  deactivated  and 
the  integrating  capacitor  values  are  recorded  to  memory  introducing  readout  noise  [2], 

The  amplifier  and  A/D  conversion  process  produce  current  values  whose  units 
are  not  in  terms  of  photo-electrons,  due  to  the  various  scaling  factors.  For  specific  am¬ 
plifier  and  reference  reverse  bias  voltage  settings,  the  system  can  be  calibrated  against 
a  known  radiometric  source  to  produce  true  photo-electron  counts.  The  gain  variation 
phenomenon  significantly  complicates  this  calibration  procedure.  Fortunately,  having 
this  exact  relationship  is  not  necessary  for  ranging. 


25 


3.6  Model  Implications 

The  first  major  implication  of  the  proposed  system  model  is  nonlinearity.  If  one 
record  is  taken  under  ambient  lighting  conditions,  a  certain  photonic  load  is  induced 
on  the  system,  and  the  reverse  bias  voltage  varies.  If  a  second  record  is  taken  at 
a  brighter  ambient  condition,  a  greater  photonic  load  is  incident,  and  the  reverse 
bias  voltage  varies  in  a  different  fashion.  Since  Miller’s  equation  is  nonlinear  with 
respect  to  reverse  bias  voltage,  the  difference  between  the  two  records  is  not  the 
linear  difference  of  the  incident  light  levels.  Even  dark  current  subtraction  must  be 
treated  with  caution,  because  the  dark  current  is  also  a  function  of  the  APD  gain. 
However,  with  incident  photons  only  due  to  self  emission  of  the  optics  and  a  lens  cap, 
it  is  the  most  linear  signal  which  can  be  measured  under  held  conditions. 

For  any  ranging  situation  of  interest  in  a  military  application,  there  are  always 
at  least  two  inputs.  The  ambient  light  from  the  scene  sets  the  reverse  bias  voltage, 
and  thus  the  gain,  in  motion  prior  to  any  returning  laser  pulse’s  arrival.  Even  if  the 
gain  behavior  during  the  pulse’s  arrival  is  completely  dominated  by  the  detector’s 
response  to  the  pulse,  the  initial  gain  value  is  not  constant  in  range.  Therefore,  the 
system’s  response  to  the  laser  pulse  is  not  shift-invariant. 

The  resulting  nonlinear  time-variant  ranging  response  complicates  both  data 
processing  and  ranging.  If  the  dark  response  is  taken  as  a  linear  value,  it  can  be 
subtracted  from  a  pulse  illuminated  response,  in  order  to  remove  additive  fixed  pattern 
noise,  but  the  ambient  response  cannot  be  subtracted.  The  effects  of  the  system  wide 
gain  fluctuation  due  to  optical  loading  must  also  be  removed  from  pulse  illuminated 
records,  to  produce  a  linear  time-invariant  signal  for  range  estimation. 

3.7  Residual  System  Identification 

A  mathematical  model  is  available  for  the  nonlinear  average  APD  gain,  and 
a  second-order  closed-loop  PID  model  is  proposed  for  the  voltage  regulator  plant. 
However,  no  information  is  available  regarding  the  actual  PID  parameters  of  the 
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voltage  regulator  or  the  presumed  RC  network  resulting  from  the  biasing  arrangement 
of  the  ROIC.  Additional  effects  from  the  analog-to-digital  conversion  are  intrinsically 
recorded  in  the  system’s  output,  which  must  also  be  removed  prior  to  considering 
the  system’s  response.  Since  the  system  cannot  be  disassembled  to  directly  measure 
component  responses,  a  system  identification  problem  remains. 

An  iterative  approach  is  used  to  sequentially  remove  the  effects  of  the  ROIC, 
time  variance,  and  nonlinearity  to  examine  the  closed  loop  relationship  between  the 
current  demand  on  the  voltage  regulator  and  the  reverse  bias  voltage  provided.  An 
impulse  response  for  the  voltage  regulator  is  the  desired  result,  if  the  current  demand 
and  bias  voltage  are  known  [13].  Once  the  response  of  the  voltage  regulator  is  known, 
it  can  be  coupled  with  Miller’s  equation  and  an  iterative  process  to  predict  the  sys¬ 
tem’s  response  to  a  known  photonic  input  which  can  vary  in  time  on  a  pixel-by-pixel 
basis. 
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IV.  Model  Development  and  Intermediate  Findings 

4-1  Introduction 

This  chapter  describes  the  process  of  exposing  and  identifying  the  underlying 
voltage  regulator  behavior,  which  is  responsible  for  the  APD  gain  variation  and  time- 
variant  system  behavior.  Intermediate  system  response  findings  are  examined,  and 
data  analysis  process  improvements  are  suggested  where  appropriate.  The  experiment 
is  presented  in  its  signal  flow  order,  in  this  case  from  the  noisy  recorded  output  back 
to  the  photon  detection  process. 

4-2  Random  Error  Suppression 

The  device  manufacturer,  ASC,  suggests  that  the  dominant  random  noise  source 
in  the  system  is  read-out  noise,  which  is  typically  GWAN  [2].  This  assertion  is  sup¬ 
ported  by  the  537-Mhz  sampling  rate  and  associated  analog-to-digital  conversion  pro¬ 
cess  to  reach  200-Hz  image  sampling  across  the  array  [1].  This  random  noise  distri¬ 
bution  can  be  observed  without  consideration  of  the  proposed  nonlinear  time- variant 
model. 

The  probability  distribution  shown  in  Figure  4.1  is  derived  from  the  presumption 
that  the  mean  of  one  hundred  records  of  the  same  range  slice  represents  the  true 
image,  and  deviations  from  this  value  at  each  pixel  represent  random  noise.  Since 
this  frame  is  an  actual  laser  pulse  return,  it  is  expected  to  have  a  contribution  from 
the  Poisson  photonic  detection  process  and  the  GWAN  digitization  and  read-out 
process  [1].  While  this  distribution  is  not  a  perfect  fit  to  a  Gaussian,  the  Gaussian 
form  appears  to  dominate  the  noise  behavior.  Were  this  not  the  case,  the  symmetry  of 
the  distribution  would  be  further  distorted,  unless  the  Poisson  lambda  value  exactly 
equalled  the  mean  GWAN  value. 

Multiple  frame  averaging  is  used  to  suppress  the  random  error,  from  all  sources, 
in  the  data  collected  for  this  work.  While  other  filtering  methods  are  necessary  for 
military  applications  when  the  target  cannot  be  imaged  multiple  times,  these  are  not 
pertinent  to  the  development  of  the  model.  The  signal  to  noise  ratio  can  be  improved 


Figure  4.1:  Illuminated  scene  noise  distribution 

by  a  factor  of  ten  for  signals  with  GWAN  by  simply  averaging  one  hundred  frames. 
This  amount  of  averaging  appears  to  suppress  random  error  effects  sufficiently  to 
proceed  with  systematic  error  suppression. 

4-3  ROIC  Systematic  Error  Suppression 

The  systematic  error  contribution  of  the  ROIC  is  generally  modeled  as  a  mul¬ 
tiplicative  gain  variation  between  amplifiers  and  an  additive  offset  for  each  amplifier. 
A  multiplicative  and  additive  response  is  also  present  due  to  the  individual  APDs  in 
the  arrays.  Because  APD  responsivity  variation  is  under  investigation,  amplifier  gain 
variations  cannot  be  evaluated  at  this  stage  of  the  model.  However,  the  additive  offset 
terms  can  be  examined  when  the  APD  multiplicative  term  is  acting  on  a  near-zero 
value. 

In  order  to  introduce  a  minimum  of  photonic  input,  a  series  of  frames  was  taken 
with  the  lens  cap  on  the  flash  LADAR’s  optical  system.  Even  in  this  dark  condition, 
some  electrons  will  be  released  by  the  APD  due  to  photons  from  self-emission  of  the 
optical  components  and  casing,  as  well  as  thermal  effects  within  the  APD  itself  [3]. 


29 


As  previously  indicated,  the  mean  individual  APD  dark  current  was  14.4-nA  prior 
to  read-out  amplification  for  a  developmental  version  of  this  device  [2],  The  dark 
response  current  demand  is  as  low  as  it  can  be  made  for  the  system,  and  the  smallest 
possible  gain  variation  due  to  optical  loading  is  experienced.  Therefore,  the  results 
will  be  attributed  to  the  ROIC  as  an  additive  component. 

The  additive  components  from  the  dark  response  represent  fixed  pattern  noise 
in  the  data  cube,  the  values  of  which  vary  both  by  pixel  and  by  range  sample.  A  series 
of  dark  response  measurements  at  increasing  range  delays  was  taken  and  the  averaged 
results  for  one  pixel  with  error  bounds  are  shown  in  Figure  4.2.  The  pattern  of  the 
recorded  value  is  clearly  repeated,  but  its  amplitude  increases  with  increasing  range. 
During  the  recording  time,  the  incident  photons  and  dark  current  from  the  APDs 
cannot  be  fluctuating  in  a  consistent  enough  manner  to  demonstrate  this  behavior 
after  averaging.  If  this  was  the  result  of  thermal  build-up,  the  responsivity  would 
have  decreased  as  the  system  warmed  up,  rather  than  increased.  Further  justifica¬ 
tion  for  assigning  the  additive  offset  values  to  the  ROIC  can  be  found  by  examining 
the  variance  around  the  mean  values.  Once  again,  a  relative  frequency  probability 
distribution  indicative  of  GWAN  is  observed  in  Figure  4.3,  with  even  less  symmetry 
distortion  than  the  illuminated  case  shown  in  Figure  4.1. 

For  the  dark  response  case,  the  APD  gain  has  been  assumed  to  be  as  close  to 
constant  valued  and  LTI  as  it  can  be  in  operation  in  this  system.  The  additive  fixed 
pattern  noise  attributed  to  ROIC  is  also  assumed  to  be  the  product  of  an  LTI  process. 
Therefore,  any  illuminated  data  cube  represents  a  sum  of  the  fixed  pattern  result  and 
the  output  of  the  non-LTI  system  response  under  optical  loading.  If  this  is  so,  it  is 
mathematically  valid  to  subtract  fixed  pattern  noise,  which  never  entered  a  nonlinear 
process,  from  an  illuminated  result. 

It  is  important  to  note  that  it  is  not  strictly  appropriate  to  subtract  an  illumi¬ 
nated  data  cube  from  another  illuminated  cube,  such  as  a  record  of  the  same  day-lit 
scene  with  and  without  a  laser  pulse  present.  When  the  pulse  is  active,  the  laser  pulse 
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Figure  4.2:  Dark  response  noise  structure  and  standard  deviation 


x  10  3 


Figure  4.3:  Dark  response  noise  distribution 
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and  ambient  photons  are  linearly  added  before  entering  the  nonlinear  time-varying 
process,  then  the  fixed  pattern  noise  is  added  to  the  result.  If  an  ambient  data  cube  is 
subtracted  from  a  pulse  illuminated  data  cube,  the  additive  pattern  noise  is  removed, 
but  the  remaining  signal  is  deformed  by  subtracting  the  final  results  of  two  nonlinear 
time-varying  processes.  Therefore,  only  dark  responses,  or  ambient  responses  ap¬ 
proaching  the  dark  response,  should  be  directly  subtracted  from  illuminated  data  to 
remove  additive  fixed  pattern  noise. 

4-4  Time  Variance 

The  time  variance  of  the  system  is  the  next  obstacle  to  be  overcome,  in  order 
to  examine  the  voltage  regulator  response.  The  proposed  model  indicates  a  system 
which  is  time  variant,  in  that  each  pixel  contributes  to  the  system’s  current  demand 
according  to  when  photons  are  incident  upon  that  pixel.  Likewise,  the  scale  of  each 
pixel’s  contribution  to  the  total  current  demand  depends  upon  the  gain  factor,  which 
varies  in  time  according  to  the  total  current  demand  and  system  response. 

Since  the  time  variant  response  arises  from  pixels  demanding  photocurrent  at 
different  times,  the  time  variance  can  be  suppressed  by  exposing  the  entire  array  to 
an  input  with  uniform  time  varying  behavior.  A  step  input  across  the  entire  APD 
array  is  a  simple  and  direct  method  of  achieving  this  result.  If  the  APD  array  is 
exposed  to  a  constant  irradiance,  but  no  reverse  bias  voltage  is  applied,  no  photon 
driven  output  response  is  expected.  As  soon  as  bias  voltage  is  applied,  the  output  will 
begin  responding  to  the  irradiance  on  the  detector.  Thus,  the  application  of  reverse 
bias  voltage  effectively  generates  an  approximate  step  input  response  from  a  constant 
irradiance. 

A  halogen  lamp  with  a  hemispheric  reflector  was  used  to  provide  an  approx¬ 
imately  constant  photon  exposure  on  each  pixel.  Radiometry  indicates  that  there 
will  be  variations  in  irradiance  between  pixels  due  to  the  physical  orientation  of  the 
lamp  and  the  flat  array  [3].  ffowever,  with  the  time  variance  suppressed,  the  average 
response  of  the  APD  array  to  an  average  irradiance  can  be  determined  as  if  it  were 
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a  single  time  invariant  detector.  This  result  allows  the  system  response  to  be  further 
evaluated. 

4-5  Time  Concatenation 

Figure  4.4  indicates  that  the  APD  array  is  biased  at  the  initial  activation  for  each 
frame,  rather  than  just  prior  to  recording.  Therefore,  concatenating  successive  records 
can  be  exploited  as  a  method  to  record  a  longer  response  than  a  single  frame,  and 
provide  a  much  more  accurate  overall  system  response  measurement.  Concatenation 
is  shown  not  only  to  be  viable,  but  to  indicate  another  source  of  systematic  error 
which  otherwise  would  not  have  been  detectable. 

A  series  of  fifty-four  dark  responses  at  successively  increasing  range  delays  and 
fifty-four  corresponding  illuminated  responses  were  recorded.  One  hundred  frames 
were  taken  for  each  record  and  averaged  to  improve  the  signal-to-noise  ratio.  Each 
dark  response  record  was  subtracted  from  its  corresponding  illuminated  record  and 
the  results  re-ordered  to  reflect  their  proper  temporal  order. 

Each  record  has  an  associated  start  delay  distance  which  is  recorded  down  to 
one-quarter  foot  intervals.  For  alignment,  each  record  was  up-sampled  by  a  factor  of 
four  and  shifted  by  integer  increments  to  its  indicated  start  location.  The  results  are 
shown  in  Figure  4.5. 

As  expected,  each  pixel’s  records  did  not  overlay  perfectly,  as  shown  for  three 
records  in  Figure  4.6.  The  beginning  of  each  successive  record  appears  to  be  attenu¬ 
ated  compared  to  the  previous  record  at  the  same  range,  but  agrees  well  later  in  its 
record.  However,  the  overall  form  of  the  overlaid  records  indicates  that  concatena¬ 
tion  is  possible,  if  a  peak  value  envelope  is  taken  to  represent  the  signal’s  true  form. 
Beginning  with  the  second  record  for  each  pixel,  the  maximum  value  of  each  record 
was  retained  in  its  area  of  overlap  with  the  following  record. 

The  averaged  concatenated  record  from  1089  pixels  near  the  center  of  the  array 
was  taken  to  represent  the  average  system  response  to  a  step  input.  This  average 
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Figure  4.4:  Halogen  lamp  illuminated  responses  minus  dark  respones  for  ranges  a. 
0-ft  b.  9. 25- ft  c.  18.5-ft  d.  27.75-ft 


Figure  4.5:  Overlay  of  illuminated  conditioned  responses  for  fifty-four  ranges. 
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Figure  4.6:  Overlay  of  three  conditioned  illuminated  response  records 

result,  seen  in  Figure  4.7,  shows  improved  smoothness  when  compared  to  individual 
records,  due  to  the  resulting  attenuation  of  any  discontinuities  in  the  overlap  region 
specific  to  one  pixel’s  set  of  records.  After  further  interpretation,  this  signal  was  used 
to  determine  the  voltage  regulator  response. 

4-6  Overlap  Deviation 

The  consistency  of  the  overlap  error  between  records  prompted  further  exam¬ 
ination  of  the  effect.  In  the  absence  of  concatenation,  no  indication  can  be  had  as 
to  what  the  expected  value  at  a  certain  range  sample  should  be.  However  with  the 
signals  concatenated,  it  is  hypothesized  that  the  trailing  portion  of  each  record  rep¬ 
resents  the  true  value.  Since  additive  fixed  pattern  noise  from  the  ROIC  has  already 
been  removed,  the  deviations  of  the  leading  portion  of  the  records  are  hypothesized  to 
be  the  results  of  a  multiplicative  process.  Based  on  this  hypothesis,  Figure  4.8  shows 
an  average  multiplicative  correction  factor  for  the  overlapping  region,  as  well  as  the 
standard  deviation  associated  with  this  average  at  each  point. 
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Figure  4.7:  Concatenated  record  maximum  value  envelope 

In  light  of  the  dark  current’s  increasing  baseline  value  over  time,  an  analysis 
of  the  temporal  behavior  of  the  overlap  deviation  was  also  conducted.  No  trend  was 
identified  regarding  the  behavior  earlier  or  later  in  the  sequence,  as  it  was  for  the 
recorded  dark  current  values.  Rather,  the  overlap  deviations  of  those  range  samples 
which  occurred  prior  to  the  expected  start  time  of  the  step-like  input  were  found  to 
have  a  higher  variance  and  behavior  less  consistent  with  the  averaged  values  than 
range  segments  which  occurred  later. 

Due  to  the  consistency  of  this  factor  through  1089  pixels  and  53  overlapping 
regions,  it  is  recommended  that  this  factor  be  applied  to  recorded  data  even  in  the 
absence  of  record  concatenation.  If  it  is  not  applied,  the  beginning  of  a  record  will 
be  attenuated  and  any  pulse  return  shape  modified  by  the  inverse  of  the  correction 
factor.  This  type  of  deviation  will  cause  most  sub-sampling  time  increment  range  de¬ 
termination  algorithm  estimates  to  indicate  an  object  more  distant  from  the  LADAR 
than  its  actual  range. 


36 


Figure  4.8:  Correction  factor  for  overlap  with  standard  deviations 

4-7  Time  Domain  Data  Interpretation 

Thus  far,  an  output  response  to  an  approximate  photonic  step  input  has  been 
conditioned  and  concatenated.  This  result  is  re-interpreted  to  allow  the  nonlinear 
process  to  be  removed  and  find  a  linear  system  relationship  between  the  current 
demand  and  the  reverse  bias  voltage.  The  following  process  was  used  to  determine 
the  reverse  bias  voltage  and  current  demand  during  the  record. 

The  output  record  represents  the  charge  stored  in  the  integrating  capacitors 
after  analog-to-digital  conversion,  and  therefore  goes  as  the  current  demand  placed 
on  the  voltage  regulator.  This  result  is  best  kept  in  its  unmodified  form,  because 
the  analog-to-digital  converter  introduces  scaling  and  includes  lower  and  upper  bound 
limits  which  can  significantly  effect  how  a  given  capacitor  charge  is  interpreted.  It  also 
means  that  this  result  must  be  re-derived  if  the  analog-to-digital  conversion  settings 
are  changed. 

The  photon  input  demand  is  presumed  constant  throughout  the  record,  since 
the  detector  was  illuminated  by  a  lamp  rather  than  a  pulse.  In  accordance  with  the 
proposed  system  model,  all  of  the  system  feedback  is  presumed  to  enter  the  recorded 
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output  via  gain  variation.  Therefore,  if  properly  scaled,  the  result  shown  in  Figure 
4.7  represents  the  gain  variation  in  time,  when  the  system  is  exposed  to  a  photonic 
step  input. 

A  digital  count  value  must  be  chosen  which  represents  the  system  operating  at 
its  designed  gain  value,  and  several  interpretations  are  possible,  depending  upon  the 
voltage  regulators  startup  behavior.  It  is  possible  that  the  gain  oscillates  around  the 
desired  value  near  the  end  of  the  record  shown  in  Figure  4.7.  It  is  also  possible  that 
the  rise  time  of  the  voltage  was  fast  enough  that  the  slight  knee  in  the  first  oscillation 
at  sample  two  actually  represents  the  desired  gain  and  that  the  system  could  not 
deliver  enough  current  to  maintain  the  desired  value  in  the  steady  state.  Having  no 
compelling  evidence  at  this  point  for  either  case,  both  are  retained  for  evaluation. 

If  the  reverse  breakdown  voltage  for  the  APDs  is  consistent  and  known,  equation 
(3.1)  can  be  solved  for  the  reverse  bias  voltage.  While  Sensors  Unlimited  indicated 
very  consistent  reverse  breakdown  voltages  during  the  development  of  the  APD  array, 
a  definitive  record  of  the  reverse  breakdown  voltage  of  the  final  device  is  not  available 
[2] .  Lacking  this  information,  the  reverse  breakdown  voltage  listed  for  a  developmental 
device  was  used  as  an  estimated  value. 

Both  possible  reverse  bias  voltages  and  both  postulated  desired  gain  scaling 
points  were  examined.  It  was  found  that  using  the  software  set  value  with  the  desired 
gain  point  set  at  the  knee  resulted  in  a  modeled  reverse  bias  voltage  with  an  imaginary 
component,  when  the  gain  fell  below  one.  With  the  reference  voltage  set  one  volt  below 
reverse  breakdown,  the  reverse  bias  voltage  model  yielded  imaginary  values  when  the 
reverse  bias  voltage  became  greater  than  the  reverse  breakdown  voltage.  Therefore, 
only  two  of  the  permutations  were  evaluated.  In  the  first  case,  the  reverse  bias  voltage 
was  set  one  volt  below  the  reverse  breakdown  voltage,  with  the  gain  point  set  at  the 
knee.  In  the  second  case,  the  reverse  bias  voltage  was  set  at  the  software  recorded 
value  with  the  gain  point  set  at  the  mean  of  the  final  two  hundred  samples  considered. 
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Using  the  developmental  reverse  breakdown  voltage,  equation  (3.1)  is  solved  for 
the  reverse  bias  voltage,  as  shown  in  equation  (4.1).  To  strictly  indicate  the  control 
relationship  between  the  reference,  Vref,  and  the  actual  reverse  bias  voltages,  the 
reverse  bias  voltage  is  redefined  as  the  reference  voltage  minus  a  voltage  variation, 
Vd,  in  equation  (4.2).  Each  of  the  postulated  gain  variation  records  is  then  used  as 
an  input  to  equation  (4.2)  yielding  the  reverse  bias  voltage  deviation  which  must  have 
given  rise  to  the  recorded  gain  variation. 

After  scaling  the  data  to  a  gain  profile,  calculating  the  voltage  deviation  re¬ 
quired  to  produce  that  gain,  and  assigning  the  current  demand  to  be  proportional  to 
the  digital  count,  two  signals  for  each  case  are  available  for  analysis.  The  resulting 
current  demand  and  the  reverse  bias  voltage  deviations  under  consideration  are  shown 
in  Figure  4.9.  These  results  represent  the  input  and  output  of  the  voltage  regulator, 
presumably  a  linear  time  invariant  system.  Any  transfer  function  or  difference  equa¬ 
tion  resulting  from  theses  values  must  be  used  with  a  certain  caution,  because  the 
input  signal  has  limited  frequency  content.  For  instance,  the  high  frequency  content 
of  an  impulse  or  step  response  cannot  be  recreated  accurately  from  these  signals.  For¬ 
tunately,  the  signals  of  interest  are  the  convolution  of  a  smoothly  varying  pulse  shape 
with  a  target  shape  which  will  tend  to  elongate  the  return.  With  this  caveat,  the 
signals  were  used  as  input  and  output  conditions  for  a  system  identification  process. 


VD  =  Vref  -  V, 


'BD 


4-8  Linear  Constant- Coefficient  Difference  Equation 

Since  the  input  and  output  signals  of  the  unknown  system  are  discrete,  and 
a  significant  time  history  is  available,  linear  constant-coefficient  difference  equations 
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Figure  4.9:  a.  Case  1  current  demand  b.  Case  1  voltage  deviation  c.  Case  2  current 
demand  d.  Case  2  voltage  deviation 


(LCCDE)  can  be  used  for  system  identification.  The  difference  equations  take  the 
form  of  equation  (4.3)  where  x  terms  are  input  values  and  y  terms  are  output  values  at 
successively  delayed  indexes.  This  relationship  is  solved  for  equation  (4.4)  to  predict 
the  next  output  value  from  a  series  of  known  past  input  and  output  values  [11]. 


x{n)  +  b\x{n  —  1)  +  ...  +  bqx(n  —  q)  —  y(n)  +  a0y(n  —  1)  +  ...  +  ary(n  —  r)  (4.3) 


y(n)  =  x(n)  +  b\x{n  —  1)  +  ...  +  bqx(n  —  q)  —  ao  y(n  —  1)  —  ...  —  ary(n  —  r )  (4.4) 

An  entire  set  of  coefficient  values  is  found  by  solving  a  system  of  LCCDEs 
relating  past  known  input  and  output  values  to  predict  the  current  output  value 
using  a  fixed  number  of  input  and  output  coefficients.  Initial  conditions,  for  the 
system  of  equations,  are  taken  as  zero  for  the  current  demand  and  reverse  bias  voltage 
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deviation  values.  The  number  of  coefficients  indicates  the  overall  order  of  the  system. 
Algebraically,  this  would  lead  to  a  great  deal  of  manipulation  for  a  large  number  of 
coefficients. 

Linear  algebra  can  be  used  to  simplify  the  solution  by  preparing  an  appended 
matrix  of  input  terms,  output  terms  and  an  identity  matrix  and  then  finding  the 
reduced  row  echelon  form  of  the  entire  matrix  [11],  Equation  (4.5)  shows  an  example 
with  three  input  coefficients  in  x  and  two  output  coefficients  in  y.  The  sixth  column 
contains  the  current  output  values,  in  this  example. 
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After  solving  for  the  reduced  row  echelon  form,  the  coefficient  values  for  each 
position  are  found  in  the  place  of  the  previous  current  output  values.  In  this  example, 
these  values  appear  in  the  sixth  column.  The  appended  identity  matrix  now  contains 
the  inverse  matrix  of  the  input  and  output  values  [11].  This  result  is  shown  in  equation 
(4.6). 
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4-9  LCCDE  Solution 


Since  the  system  order  is  unknown  and  the  solution  time  for  a  given  matrix 
is  short,  an  error  optimization  method  was  used  to  determine  how  many  coefficients 
should  be  used.  The  matrix  solution  was  found  for  one  to  two-hundred  input  co¬ 
efficients  and  one  to  two-hundred  output  coefficients,  then  the  known  input  values 
were  processed  recursively  by  the  derived  difference  equation.  The  error  between  the 
recorded  output  voltage  signal  and  the  recursively  derived  output  signal  was  found 
and  recorded.  Finally,  the  set  of  coefficients  which  produced  the  smallest  error  was 
selected  as  the  coefficient  set  for  the  system  model. 

The  final  output  system  model  for  Case  1  has  ninety  input  and  thirty-nine 
output  coefficients,  while  the  model  for  Case  2  has  one  hundred  and  forty-one  input 
and  fifty-five  output  coefficients.  The  results  of  reconstructing  the  reverse  bias  voltage 
deviation  records  are  shown  in  Figure  4.10. 


a  b 


Figure  4.10:  a.  Reconstructed  voltage  deviation  for  case  1  b.  Reconstructed  voltage 
deviation  for  case  2 
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V.  Model  Validation 


5.1  Forward  Process 

Based  on  the  system  response  derived  in  section  4.9,  a  direct  iterative  process 
is  used  to  predict  the  result  recorded  by  the  system  for  a  known  true  photon  arrival 
distribution  in  time,  using  the  following  assumptions.  First,  initial  conditions  of  zero 
in  voltage  variation  and  current  demand  are  assumed  for  a  set  of  samples  as  large  as 
the  largest  number  of  filter  coefficients.  Second,  the  initial  gain  is  presumed  to  be  the 
empirically  derived  value  for  the  ideally  regulated  reverse  bias  voltage. 

The  current  demand  at  the  time  index,  k,  is  determined  by  multiplying  the  sum 
of  the  photon  arrival  count,  P(q,r,k),  for  each  pixel  at  the  time  index  by  the  gain, 
G(k),  at  the  time  index. 


128  128 


(5.1) 


Then,  the  voltage  variation  for  the  time  index  which  arises  from  all  previous 
voltage  variations  and  current  demands  is  recursively  calculated  using  equation  (5.2). 


VD(k )  =  ID{k)  +  biID{k-l)  +  ...  +  bzID{k- z)-a1VD(k-l)~  ...  -  azVD(k  -  z)  (5.2) 


The  resulting  voltage  variation  at  the  time  index  is  then  used  in  equation  (3.1) 
to  determine  the  gain  for  the  next  time  index.  This  process  is  repeated  until  the 
known  photon  arrival  values  are  exhausted,  and  a  gain  profile  has  been  calculated  for 
the  entire  record  length.  Given  that  only  a  forward  process  is  proposed,  a  limited  test 
case  is  required  for  which  the  true  photon  arrival  distribution  is  known. 

5. 2  Forward  Simulation 

If  the  system  model  is  valid  over  its  range  of  assumptions,  the  best  measure  of 
its  performance  is  its  ability  to  recreate  the  gain  variation  experienced  by  a  simple  test 
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case.  Unfortunately,  the  test  cases  available  for  comparison  were  limited  to  data  taken 
on  previous  collections.  Since  no  cases  were  available  for  a  single  surface  occupying  the 
entire  FOV,  a  compromise  case  was  selected.  The  derived  noise  reduction  procedures 
were  applied  to  the  recorded  data  and  the  RGB  estimation  method  was  used  to 
determine  range  values  using  the  original  and  modified  system  models. 

The  test  scenario  consists  of  a  two-surface  target  boards,  the  ground,  the  am¬ 
bient  lighting  conditions,  and  some  foliage  near  the  end  of  the  range  record.  While 
not  ideal,  this  data  set  was  taken  at  night,  minimizing  the  gain  variations  due  to  the 
system’s  initial  step  response  when  exposed  to  ambient  light.  The  return  from  the 
ground  and  vegetation  causes  some  deviation  from  the  modeled  case,  and  the  laser 
pulse  shape  is  not  precisely  recorded.  An  additional  obstacle  is  the  lack  of  a  record 
of  the  target  surface  separation  or  true  target  range  for  this  data  set. 

Following  the  procedure  for  additive  fixed  pattern  noise  removal  developed  in 
section  4.5,  a  dark  response  measurement  for  the  range  delay  of  the  test  data  set  was 
taken.  Both  the  dark  response  and  recorded  data  set  consisted  of  one-hundred  frames. 
Therefore,  random  noise  suppression  on  par  with  the  model  development  values  was 
achieved  through  frame  averaging.  The  averaged  result  of  the  dark  response  was 
subtracted  from  the  averaged  illuminated  frame  to  generate  the  final  frame  for  the 
ranging  algorithm. 

The  RGB  estimation  theory  method  was  selected  as  the  ranging  algorithm  for 
the  test.  This  selection  was  due  to  its  simplicity  of  implementation  and  relevance 
to  ongoing  research.  The  algorithm  was  run  on  the  recorded  data  set  in  its  original 
form  with  a  resolution  of  1/30  of  a  range  sample,  or  approximately  2-cm,  and  on  the 
recorded  data  modified  by  a  varying  gain  profile  at  the  same  resolution.  To  derive  a 
gain  profile,  each  pair  of  postulated  inverted  parabola  amplitudes  and  locations  was 
used  as  a  true  photonic  input  into  the  non-LTI  system  model.  Thus,  the  gain  profile 
was  different  in  each  iteration,  depending  upon  the  hypothesized  ranges  of  the  target 
surfaces. 
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The  photonic  scaling  for  the  input  parabolas  was  somewhat  arbitrary,  being 
based  on  previously  calculated  radiometric  values  for  the  target.  For  the  purposes  of 
the  model,  as  long  as  the  photonic  pulse  amplitude  is  the  correct  factor  greater  than 
the  bias  photonic  level,  the  scaling  factor  in  the  ranging  algorithm  will  rescale  the 
results  to  match  the  digital  count  record.  The  photonic  bias  level  was  taken  from  the 
background  recorded  mean  and  the  photonic  pulse  scaled  accordingly.  Since  the  front 
and  back  surface  areas  of  the  test  target  exposed  to  the  LADAR  were  of  approximately 
equal  area,  the  total  photon  input  was  modeled  as  the  sum  of  two  pulses.  This  total 
photon  input  was  run  through  the  forward  process  described  in  section  5.1,  and  the 
results  used  to  calculate  the  MSE  for  each  postulated  location. 

The  algorithm  was  run  for  both  the  Case  1  and  Case  2  coefficient  sets.  It  was 
found  that  using  Case  2  coefficients  produced  an  initial  voltage  variation  due  to  the 
Erst  surface  pulse  which  generated  imaginary  gain  profile  values.  Therefore,  Case  2 
was  dropped  from  consideration.  The  results  for  Case  1  are  presented  in  Figure  5.1. 
It  is  important  to  note  that,  while  Case  2  was  dropped,  the  Case  1  coefficients  have 
not  been  proven  to  be  absolutely  correct.  Rather,  the  Case  1  coefficients  are  shown 
to  be  the  best  model  for  the  actual  system  response,  of  those  under  consideration. 

A  comparison  of  the  pulse  center  locations,  for  Case  1,  gives  the  following  results. 
The  unmodified  RGB  estimation  theory  method  indicated  a  first  surface  location  of 
range  sample  6.36  and  a  second  surface  location  of  9.06,  based  on  a  minimum  MSE 
of  5.747  *  103-digital  counts  and  1.755  *  104-digital  counts,  respectively.  The  RGB 
estimation  theory  algorithm  using  the  proposed  system  gain  model  in  lieu  of  the 
original  bias  estimation  method  yields  a  first  surface  location  of  range  sample  6.46 
and  a  second  surface  location  of  9.16,  with  minimum  MSEs  of  3.922  *  103-digital 
counts  and  5.955  *  103-digital  counts.  The  gain  model  system  results  in  a  factor  of 
1.465  improvement  in  the  minimum  MSE  for  the  first  surface,  and  a  factor  of  2.946 
improvement  in  the  minimum  MSE  for  the  second  surface. 
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Range  Sample 


Range  Sample 


Figure  5.1:  a.  Original  model  and  data  for  first  surface  b.  Original  model  and  data 
for  second  snface  c.  Modified  model  and  data  for  first  surface  d.  Modified  model  and 
data  for  second  surface 


The  results  from  the  modified  estimation  theory  method  predict  surface  loca¬ 
tions  more  distant  from  the  LADAR  than  the  unmodified  method.  This  is  in  line  with 
the  expectation  that  falling  gain  as  a  pulse  arrives  will  tend  to  skew  the  pulse  center 
away  from  the  detector.  Even  handicapped  by  the  inability  to  postulate  the  gain 
variation  due  to  the  return  from  the  ground  prior  to  the  target,  the  MSE  reductions 
indicate  a  significant  model  improvement.  This  is  particularly  true  for  the  second 
surface,  for  which  the  gain  behavior  can  be  reasonably  expected  to  be  dominated  by 
the  first  surface  return’s  effects.  Pulse  reshaping  is  also  evident  in  the  closer  fit  of  the 
corrected  pulse  shape  to  the  measured  data,  when  compared  to  the  uncorrected  pulse 
shape. 

This  experiment  cannot  show,  due  to  the  lack  of  an  independent  range  mea¬ 
surement,  a  definitive  improvement  in  the  range  estimate  for  the  test  case.  However, 
it  does  show  a  significant  reduction  in  MSE  between  the  recorded  data,  and  the 
modeled  pulse  which  determined  the  range  estimate.  At  a  minimum,  the  resulting 
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improvement  of  the  fit  between  the  modeled  and  recorded  data  indicates  the  validity 
of  incorporating  a  varying  gain  component  in  the  system  model. 
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VI.  Conclusion  and  Recommendations 


6.1  Primary  Areas  of  Contribution 

The  primary  contribution  of  this  work  is  to  explain  and  model  the  anomylous 
behavior  found  in  data  from  the  AFLR/SNJM  FLASH  LADAR  as  a  non-LTI  gain 
variation  process.  It  was  determined  that  background  subtraction  in  illuminated 
scenes  must  be  used  cautiously,  due  to  non-LTI  nature  of  the  system.  The  use  of  dark 
current  subtraction  to  remove  the  additive  offset  component  of  each  pixel  and  A/D 
converter  has  been  proposed  as  a  means  to  reduce  additive  fixed  pattern  noise  while 
introducing  minimal  error.  Finally,  a  responsivity  variation,  which  is  attributed  to 
the  read-out  amplifier,  is  identified  and  characterized. 

6.2  Significance  of  Contribution 

The  significance  of  the  forward  model  cannot  be  understated  in  regards  to  sub¬ 
sample  accuracy  range  determination  algorithms,  which  are  derived  for  LTI  signals. 
The  gain  variation  introduces  dependence  between  all  of  the  pixels  in  the  array,  which 
was  neither  expected  nor  modeled  previously.  In  the  face  of  pixel  interdependency  due 
to  gain  variation  and  the  optical  system’s  PSF,  the  achievable  accuracy  and  resolution 
of  such  algorithms  can  only  approach  the  analytically  derived  values  in  very  limited 
cases. 

If  multiple  range  returns  are  to  be  interrogated  from  a  single  IFOV,  the  gain 
variation  must  be  removed.  Otherwise,  true  multiple  returns  may  be  obscured  by 
the  gain  variation,  and  false  secondary  returns  may  be  observed.  Removing  the  gain 
profile  will  also  alleviate  false  ranges  being  applied  to  pixels  observing  targets  outside 
the  range  gate. 

The  modification  of  the  background  subtracting  scheme  represents  a  significant 
refinement  in  the  understanding  of  the  system.  The  proposed  dark  response  method 
introduces  less  error  than  the  original  illuminated  scene  method,  but  does  not  address 
the  contributions  of  ambient  light.  To  properly  remove  the  ambient  light  contribu- 
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tions,  the  ambient  response  must  be  recorded  and  linearized  before  subtraction  from 
a  linearized  pnlse  illuminated  signal. 

The  amplitude  scaling  identified  during  record  concatenation  also  effects  ranging 
accuracy,  by  non-uniformly  rescaling  the  recorded  returns.  If  only  single  records  were 
examined,  this  deviation  would  likely  have  been  lumped  into  APD  gain  variation 
rather  than  identified  as  a  separate  source  of  error. 

6.3  Recommendations  for  Further  Work 

The  first  recommended  continuation  of  this  experiment  relates  to  the  Miller’s 
Formula  based  APD  gain  term.  The  near  steady  state  responsivity  for  a  series  of  bias 
voltages  should  be  investigated  to  refine  and  scale  Miller’s  Equation  for  this  device, 
at  fixed  amplifier  settings.  This  may  be  accomplished  by  exposing  the  detector  to  a 
series  of  known  input  irradiances  and  concatenating  a  series  of  records  at  the  longest 
delays  allowed  by  the  system.  This  will  allow  the  worst  of  the  transient  gain  effects 
to  atrophy,  and  an  average  of  the  oscillating  value  may  be  taken  as  the  desired  value. 
Making  these  measurements  will  answer  with  certainty  where  the  gain  set  point  lies 
at  the  beginning  of  the  record. 

It  may  be  possible  to  investigate  the  system’s  response  directly.  A  laser  with 
a  pulse  length  much  shorter  than  the  range  sampling  time  increment  can  be  used 
to  generate  an  impulse  like  response.  If  the  input  pulse  intensity  is  not  too  high, 
the  gain  variation  following  the  pulse  may  be  visible  even  in  dark  conditions  without 
saturation.  An  extended  record  can  be  taken  by  fixing  the  laser  pulse  timing  and 
adding  progressively  more  delay  for  recording  in  sular  mode. 

After  refining  the  gain  variation  model,  the  next  critical  step  to  be  taken  is 
to  develop  an  approach  to  remove  the  effects  of  gain  variation  from  the  data.  Two 
methods  present  themselves.  The  first  would  be  to  directly  measure  the  voltage 
regulator  value  during  recording.  If  this  value  were  sampled  by  the  ROIC,  a  refined 
Miller’s  Equation  could  be  used  to  directly  rescale  the  record.  Alternately,  an  iterative 
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method  could  be  implemented  based  on  the  forward  model  developed  in  this  work. 
However,  this  method  would  require  background  and  dark  measurements  of  the  same 
scene  prior  to  and  including  the  range  of  interest  and  will  likely  prove  impractical  for 
all  but  research  tasks. 
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Appendix  A.  Matlab®  Code 
A.l  Matlab®  Code  for  data  Process  One 

Listing  A.l:  Processlvl.m 

“/oProcesslvl  .  m  load  .  seq  data  files  from  a  fixed  directory  located 
on 

'/.line  63  and  averages  the  results  for  the  entire  array.  It  also  . 
calculates 

% basic  statistics  on  a  limited  set  controlled  by  the  prange  vector 
in 

7.  [1st  row,  last  row,  1st  col,  last  col]  format.  A  series  of  ... 
inputs  are 

70used  to  allow  for  partially  corrupted  data  sets  to  be  used. 

clear  all;  close  all;  clc; 
tic  , 

7,1.  preparation  for  data  loading  7«****************** 
disp (’ Process  set  #1’); 

N  =  1 08  ;  7,  NUMBER  of  data  files 

Frames  =  100;  7,  number  of  frames  per  file 

FnT  =  108;  7.  Total  number  of  files 

FnD=54;  7,Total  number  of  data  files.  It  is  expected  to  be  1/2  the 
total  since  a  dark 

‘/.current  value  is  needed  for  each  illuminated  file. 


gapinit  =  [0];  7«gap  in  file  positions  for  corrupted  files  -  0  for  . 
none 

gap=[gapinit  ( gapinit -FnD )  ]  ;  7.overlap  gap  vector  for  conditioning 
of  runs  which  should  not  be  processed  due  to  later  gaps 

7.  ad  j  udi  cat  e  number  of  good  files 
if  gapinit ~  =  0 ; 

FnG  =  FnT  -  length  (  gapinit  )  ;  /.Number  of  Good  files 

else 

FnG=FnT ; 

end  ; 

7,selction  of  pixel  range  under  consideration 
7,  [1st  row,  last  row,  first  col,  last  col] 
prange=[32  64  32  64] ; 

row=prange(l) : prange (2) ; 
col=prange (3) : prange (4) ; 
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%  initialize  the  pixel  value  matrix 
Pmat  =  zeros (length (row) ^length (col)  ,2)  ; 

count  =  1 ; 

for  n=l : length (row)  ; 

for  m  =  l : length ( col )  ; 

Pm at ( count  ,l)=row(n)  ; 

Pm at ( count  ,2)=col(m)  ; 
count  =  count  + 1 ; 

end  ; 

end  ; 

7«find  the  number  of  pixels  under  consideration  for  basic  ... 

statistics 
temp=size(Pmat) ; 
pixels  =  temp  (  1 )  ; 

12.  Load  data  for  overall  operations  &  segment  evaluation 

for  q  =  l:l:FnT;  '/.incrment  through  number  of  files 

if  sum (gap==q) ~=1 ;  “/.if  not  a  skip  file,  load  data 

[D,H, stop, mark]=flash_read_seqv2( ’d:\seal\d at a  10-2-06\’ ,[... 

’ seall00206run ’  num2str(q)  ’ .seq’]) ; 
disp (num2str (q) ) ; 

“/.initialize  matrixes  for  this  file 
Dout (q)  . data  =  zeros (128 , 128 ,20)  ; 

Dout(q)  . varMat  =  zeros (pixels  ,20)  ; 

Dout(q)  . meanMat  =  zeros (pixels  ,20)  ; 
statdata  =  zeros (pixels  ,20, Frames)  ; 

for  r  =  l  :  1  :  Frames  ;  “/.increment  through  Frames 

Dout  (q)  .  data  =  D  (r)  .  data  +  Dout  (q)  .  data  ;  “/.add  the  data  ... 
for  this  frame  to  a  total  for  this  file 

for  n  =  l  :  1  :  pixels  ;  “/.increment  through  each  pixel  for  ... 
cons ider at  ion 

for  m  =  l:l:20;  “/.increment  thru  slices  for  each  ... 
pixel 

statdata(n,m,r)=D(r)  ,data(Pmat(n,l)  , Pmat (n , 2)  , . .  . 
m)  ;  “/.compile  statdata  matrix 

end  ; 

end  ; 

end  ; 

Dout  (q)  .  data  =  Dout  (q)  .  data  . /Frames  ;  “/.find  overall  data  mean... 
matrix 

Dout  (q)  .  stop  =  stop  ;  “/.record  stop  distance 

Dout  (q)  .  pixels  =  Pmat  ;  “/.record  map  of  pixels  for  close  ... 
examinat ion 
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Dout(q).varMat=var(  statdata,[],3  );  '/.calculate  variance 
for  each  pixel  under  consideration  down  the  frames 
Dout  (q)  .  meanMat=mean  (  statdata,3  );  '/.calculate  mean  for 
each  pixel  under  consderation  down  the  frames 


end  ; 

pack  ; 
toe  , 


end  ; 


'/.save  the  output  results 

save  processldata  Dout  N  Frames  FnT  FnG  FnD  pixels  gap  Pmat  ; 


A. 2  Mat  lab®  Code  for  data  Process  Two 

Listing  A. 2:  Process2vl.m 

'/.for  those  pixels  for  which  basic  statistics  have  been  derived  in 
processl 

'/.conduct  a  dark  current  subtraction  and  concatenate  the  values  to 
form 

"/.a  single  long  data  run 
clear  all;  close  all;  clc; 

'/.show  plots  or  not  ,  1  to  show  ,  0  to  not  show 

showplot  s  =0 ; 

'/.load  data  file 
load  processldata 

"/.initialize  a  zeros  matrix  for  each  pixel 
for  m=l : 1 : pixels ; 

p(m)  . data=zeros(FnG,20)  ; 

end  ; 

"/.load  values  for  each  pixel  into  array 
for  m=l : 1 : pixels ; 
for  n=l : 1 : FnT ; 

if  sum ( gap==n) <1 ; 

p(m) . dat a (n  , :)=squeeze(Dout(n) .data(Pmat(m,l) , Pmat (m. 

,2),:)); 

end  ; 

end  ; 

end  ; 

"/.determine  the  stop  points  of  the  files 
stoplist=zeros (1 ,FnD) ; 
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for  m= 1 : 1 : FnD ; 

if  sum ( gap ==m) <  1  ; 

’/.multiply  by  4  and  round  to  insure  integer  values  of  ... 

sampling  of 
"/.stops 

stoplist (m) = round (Dout (m)  . stop*4)  ; 

end  ; 

end  ; 

7.  initialize  a  zeros  matrix  for  each  pixel 
for  m=l : 1 : pixels ; 

q(m)  . data  =  zeros (FnD ,20)  ; 

end  ; 

7.  dark  current  subtractions 
for  m=l : 1 : pixels ; 
for  n=l : FnD ; 

if  sum ( gap==n) <1 ; 

q(m)  ,data(n,  :)=p(m)  .data(n  +  54,  : ) -p (m)  .data(n,  :)  ; 

end  ; 

end  ; 

end  ; 

clear  p ; 
pack  ; 

7.  initialize  a  zeros  matrix  for  each  re-ordered  pixel 
for  m=l : 1 : pixels ; 

r(m)  . data  =  zeros (FnD  ,19)  ; 

end  ; 

7. find  the  zero  point  &  re-order 
for  m=l : 1 : pixels ; 
for  n=l : FnD ; 

if  sum ( gap==n) <1 ; 

index=find(q(m) ,data(n, :) <=1) ; 

r(m).data(n,l:(20-index))=q(m).data(n,(  (index+1) :20)... 

)  ; 

r(m).data(n,(20-index  +  l)  :  (19)  ) =q (m)  . data (n , 1 :( index ..  . 

-D); 

end  ; 

end  ; 

end  ; 

clear  q; 
pack  ; 

7.  initialize  a  zeros  matrix  for  each  pixel  preparing  for  4x  interp 
for  m=l : 1 : pixels ; 

t (m)  . data  =  zeros (FnD ,77)  ; 

end  ; 
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%  interpolate  each  string  so  they  can  be  aligned 
for  m=l : 1 : pixels ; 
for  N=1 : FnD ; 

t(m)  .data(N,  :)  =  interpl(r(m)  .data(N,  :)  ,  [0:  .25:19]  ,  ’ pchip  1  )  ; 

end  ; 

end  ; 

'/.clear  r  ; 
pack  ; 

"/.multiply  scale  for  nanoseconds  as  a  reference 
x  = ( 1 . 862/4)  . *[stoplist(l)  : 1 :  st  opl i st ( FnD ) ]  ; 

"/.initialize  a  final  data  set  variable 
sdata  =  zeros (FnD , length (x) )  ; 

for  m=l : 1 : pixels ; 

s(m)  ,data  =  zeros(FnD,length(x))  ; 

end  ; 

“/.concatenate  the  data  sets  into  one  signal 
for  m=l : 1 : pixels ; 
for  n=l : 1 : FnD ; 

if  sum ( gap==n) <1 ; 

temp=t (m) . data (n , :) ; 

sdata(n,(stoplist(n)-5*n)+6:l:(stoplist(n)+76-5*n+6)  )... 
=  temp  ; 

"/.some  hardcoded  factors  which  should  be  consistent  ... 
with  the 

"/.overlap  of  contiguous  files  are  present  +76,  -5*  ... 

f  actor 

end  ; 

s(m)  . dat a  =  max ( sdata ( 1 : FnD ,1: length(x))  ,  []  ,1)  ; 


end  ; 

end  ; 
pack  ; 

if  showplots ==1 ; 

“/.form  out  the  images  of  the  values 

for  m=l  :  1  :  pixels  ;  “/.replace  with  pixels  for  full  run 
figure  (  1 )  ; 
subplot (2,1,1)  ; 

for  n=l:l:FnD;  “/.replace  with  FnD  for  full  run 
hold  all; 
if  sum (gap==n) <1 ; 

plot (  (stoplist(n)-5*n):l:(stoplist(n)+76-5*n),  t(... 
m)  . data (n  ,  :  )  )  ; 

end  ; 

end  ; 
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title ([’ overlays  of  pixel  =  ’  num2str(m)  ’  position  =  ’  . 

num2str (Pmat (m  ,  1) )  ’,’  num2str ( Pmat (m , 2) ) ] )  ; 

xlabel ( ’ samples ’ )  ;  ylabel ( 1  photo -electron  counts’); 
hold  off;  axis  tight;  grid  on; 
axis (  [0  1767  0  1200] )  ; 

subplot (2,1,2) ; 

plot (x , s (m) . data) ;  axis  tight;  grid  on; 

title ([ 1  concatenation  of  pixel  =  1  num2str(m)  ’  position. 

=  ’  num2str (Pmat (m  ,  1)  )  ’,’  num2str (Pmat (m  ,  2)  )  ]  )  ; 

xlabel ( ’nanoseconds  ’)  ;  ylabel (’ photo -electron  counts’); 
axis (  [0  912  0  1200] )  ; 
set(l,  ’Position’  ,  [-3  39  1280  915])  ; 

pause ; 

close  (  1 )  ; 


end  ; 


end  ; 

‘/.save  process  output 

save  process2data  s  x  t  r  pixels 


A. 3  Mat  lab®  Code  for  data  Process  Three 

Listing  A. 3:  Process3vl.m 

7»This  process  is  run  after  a  series  of  data  sets  are  concatenated 
clear  all;  close  all;  clc; 

7o determine  the  impulse  response  from  a  composite  step  response 
load  process2data 

7o  averaging  concatenated  data  sets 

d at a=zeros(l, length (s(l) .data)) ; 
for  n=l : pixels ; 

data=data+s (n) .data; 

end  ; 

data=data ./pixels ; 
save  process3data  h  data 


A. 4  Mat  lab®  Code  for  data  ABhderiveknee2 
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Listing  A. 4:  ABhderiveknee2.m 

7, control  input  simulation  including  inverse  filtering  operations 
clear  all;  close  all;  clc; 


5  load  process3data ; 
clear  h; 


70downsample  the  data  back  to  the  original  scale,  this  will  ... 
introduce 

70some  interpolated  values. 

10  dat a  =  re sample ( data , 1 , 4)  ; 

7,  set  the  input  value  as  being  the  value  after  division  by 
'/.the  initial  gain,  case  #1  in  the  main  document 
p  =5 . 48 ; 

15  Vref=XX;  7.actual  value  removed  for  public  release 
Vbd  =  XX  ;  ‘/.actual  value  removed  for  public  release 

‘/.initial  gain  calculation 
Ml = 1 / (  1-abs (Vref/Vbd) ~p  ); 

20 

‘/.setting  input  scale  to  mean  in  initial  values  &  initializing  gain... 
record 

scale=mean (data (2) ) /Ml ; 
input=[scale*ones (1 ,426)]  ; 

25  M  =  data ( 1 : 426)  . / input ( 1 : 426)  ; 

Mhold  =  M  (2  :  425)  ;  ‘/.dropping  first  artifact  value 

‘/.V  g  is  the  voltage  deviation 

Vg  =  Vref - Vbd . * (1 -1 . /Mhold)  .  ~ (1/p)  ; 

30 

Vghold=Vg ; 

7. Alist  =  [90  :  1  :  92]  ,  Blist  =  [36  :  1  :  40]  to  recreate  known  values 
Alist  =  [1 : 1 : 200]  ; 

35  Blist  =  [1:1: 200]  ; 


k  =  1 

q=i 

tic 

40 


‘/.additional  sequence  delays 


“/.increment  through  the  LCCDE  solution  in  matrix  form 


for  indexl =1 : 1 : length (Alist ) ; 

for  index2  =  l : 1 : length (Blist )  ; 
45 

Acof =Alist ( indexl ) ; 

Bcof =Blist ( index2 ) ; 


‘/.the  zeros  assignment  is  critical  for  this  method  to  work 
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'/.also  I_d  goes  as  M,  but  multiplied  back  to  the  data  ... 
values  & 

'/.assume  that  every  pixel  saw  the  average  for  the  demand 
Id= [zeros (1 , floor (max(Acof  ,Bcof)+min(Bcof  ,Acof)/2))  Mhold* .  .  . 
scale  * (128~2)  ]; 

Vg= [zeros (1 , floor (max (Acof  ,Bcof)+min(Bcof  , Acof )/2) )  Vghold  .  .  . 

] ; 

'/.initialize  the  matrix  to  be  solved  via  RREF 
Ms=zeros ( Acof +Bcof , 2* ( Acof +Bcof ) + 1 ) ; 

’/.assign  identity  matrix  values 

Ms (1  :  (Acof+Bcof)  ,  (Acof  +  Bcof+2)  :  (2* (Acof+Bcof )+l) )=eye (Acof . .  . 

+Bcof , Acof+Bcof) ; 
count  =  1 ; 

7.  assign  previous  input  &  output  values 

for  n=(k  +  max(Acof  , Bcof  )  )  :  1 :  (k  +  max(Bcof  ,Acof)+Acof+Bcof-l)  ; 

Ms ( count  ,1:  (Acof  +  Bcof+1) )  =  [ Id (n : -1: (n-Bcof  +  1) )  -l*(Vg.  .  . 

((  n-1) : -1: (n-Acof)))  Vg(n)]; 
held(count) .Ms=Ms; 
count  =  count  + 1 ; 

end  ; 

“/.solve  the  system  of  equations 
Ns=rref (Ms) ; 

“/.extract  the  coefficients 
B  =  Ns(l:Bcof  , (Acof+Bcof+1))  ; 

A=[Ns((Bcof+l)  :(Bcof  +  Acof)  ,  (Acof  +  Bcof+1))  ’  ] 

“/.regenerate  the  test  input  from  the  linear  difference  ... 
equat ions 

y2= [Vg ( 1 : Bcof )  zeros (1 , 400) ] ; 
x2= Id ; 

7.  iterative  forward  solution 
for  n=(l  +  max(Acof  , Bcof  )  )  :  (400)  ; 

y2(n)=sum(x2(n:-l:n-Bcof+l) . *B,)-sum(y2(n-l:-l:n-Acof)... 
.  *  A  1  )  ; 

end  ; 

start=length([Vg(l:Bcof)  zeros(l,floor(max(Acof ,Bcof)+min(... 
Bcof ,Acof)/2))]) ; 

“/.error  calculations  &  saving  coefficient  information 

output (q)  . error  =  mean (  ( y2 ( start  : 400) -Vg ( start  : 400) ).  ~2  ); 

output ( q) . A=A ; 

output ( q) . B=B ; 

output (q) .Acof=Acof ; 

output (q)  .Bcof=Bcof  ; 
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q=q+l ; 

disp (num2str (q) ) ; 

end  ; 


end  ; 

disp ( [ " est imates  =  "  num2str (q- 1 ) ] ) ; 

toe  , 

7.  prepared  all  saved  results  for  error  evaluation 
evals=zeros (1 , length (Alist)*length(Blist) ) ; 
f  or  n=  1  : 1 :  ( q-  1 )  ; 

evals (n)=output (n) . error ; 

end  ; 

"/.search  out  minima  using  built  in  f  unct  ions 
Eindex=f ind ( evals ==min ( eval s ) ) ; 
est=l ; 

"/.notify  the  user  if  more  than  one  minima  is  found 
if  length (Eindex >1)  ; 

disp(["found  "  num2str ( length ( Eindex ) )  "  mse  mins,  displaying  ... 

first"]); 

end  ; 

"/.reload  the  coefficients  for  the  minimum  error  set 
A=output(Eindex(est)) .A; 

B=output(Eindex(est)) .B; 

Acof=output (Eindex (est) ) . Acof ; 

Bcof=output (Eindex (est)) .Bcof ; 
mse=output (Eindex (est) ) . error ; 

7.  regenerate  the  reconstruction  results  for  display 
y2=[zeros (1 ,420)] ; 
x2= Id ; 

for  n=(l  +  max(Acof  , Bcof ) )  :  (400)  ; 

y2 (n) =sum (x2(n:-l:n-Bcof+l) .  *B  ’  ) -sum (y2(n-l:-l:n-Acof) .  *A  ’  ) ; 

end  ; 
toe  , 

disp(’N0TE:  The  test  input  has  finite  frequency  content."); 

disp("the  resulting  impulse  response  cannot  address  high  frequency... 
input  s  "  )  ; 


7.  plotting  section 
figure ( 1 ) ; 

subplot (2 , 1  ,  1 )  ;  plot (Vg ( 1 : 400) )  ;  axis  tight;  title  ("  f  iltered 
output " ) ; 
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subplot (2 , 1 , 2)  ;  plot (y2  (1 : 400) )  ;  axis  tight;  title ([’  recon  output 
of  order  A=  ’  num2str ( Acof )  ’  and  B=  ’  num2str ( Bcof ) ] ) ; 

figure (2)  ; 
subplot (2,1,1)  ; 

plot (  ( abs ( Vg ( 1 : 400 ) -y2  ( 1 : 400 ) ) )  ,  ’  r -  1  )  ;  axis  tight;  title(’abs  .. 

error  by  point1); 
subplot (2,1,2)  ; 
plot (Vg  (1  :  400) )  ;  hold  all; 
plot ( y2  (  1  :  400 ))  ;  hold  off; 

title (’ overlaid  voltage  &  reconstructed  voltage’); 

‘/.save  data  for  further  use 

save  ABvalsknee2  A  B  Acof  Bcof  Vg  y2  x2 


A. 5  Matlab®  Code  for  data  ABhderivemean2 

Listing  A. 5:  ABhderivemean2.m 

'/.control  input  simulation  including  inverse  filtering  operations 

clear  all;  close  all;  clc; 

load  process3data ; 
clear  h ; 

'/.downsample  the  data  back  to  the  original  scale,  this  will  ... 
introduce 

'/.some  interpolated  values, 
dat  a  =  re sample ( data , 1 , 4)  ; 

'/.set  the  input  value  as  being  the  value  after  division  by 
'/.the  initial  gain,  case  #1  in  the  main  document 
p=5 . 48 ; 

Vref=XX;  '/.actual  value  removed  for  public  release 
Vbd  =  XX  ;  '/.actual  value  removed  for  public  release 

Ml = 1 / (  1-abs (Vref/Vbd) ~p  ); 

'/.setting  scale  to  mean  of  final  oscillations  of  samples 
scale=mean (data (350 : 420) ) /Ml ; 
input  =  [scale*ones (1 ,426)]  ; 

M  =  data(l:426)  ./input(l:426)  ; 

Mhold  =  M (2 : 425)  ; 

“/.Vg  is  the  voltage  deviation 

Vg  =  Vref - Vbd . * (1 -1 . /Mhold)  .  ~ (1/p)  ; 

Vghold=Vg ; 
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"/, Alist  =  [137  :  1  :  141]  ,  Blist  =  [53  :  1  :  57]  to  recreate  known  values 
Alist=  [1  :  1 : 200]  ; 

Blist =  [1:1: 200] ; 
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k  =  1 

q=i 

tic 


‘/.additional  sequence  delays 


40  "/.increment  through  the  LCCDE  solution  in  matrix  form 

for  indexl =1 : 1 : length (Alist ) ; 

for  index2  =  l : 1 : length (Blist )  ; 

45  Acof =Alist ( indexl ) ; 

Bcof =Blist ( index2 ) ; 

"/.the  zeros  assignment  is  critical  for  this  method  to  work 
"/.also  I_d  goes  as  M,  but  multiplied  back  to  the  data  ... 
values  & 

50  "/.assume  that  every  pixel  saw  the  average  for  the  demand 

Id= [zeros (1 , floor (max(Acof ,Bcof)+min(Bcof ,Acof)/2))  Mhold* 
scale*(128~2)  ]; 

Vg= [zeros (1 , floor (max (Acof ,Bcof)+min(Bcof ,Acof)/2))  Vghold 

]; 

55  "/.initialize  the  matrix  to  be  solved  via  RREF 

Ms=zeros ( Acof +Bcof , 2* ( Acof +Bcof ) + 1 ) ; 

"/.assign  identity  matrix  values 

Ms(l:(Acof+Bcof)  ,(Acof  +  Bcof+2)  :(2*(Acof+Bcof)+l))=eye(Acof 
+Bcof ,Acof+Bcof) ; 
count  =  1 ; 

60 

"/.assign  previous  input  &  output  values 

for  n=(k  +  max(Acof  , Bcof ) )  :  1 :  (k  +  max(Bcof  ,Acof)+Acof+Bcof-l)  ; 
Ms ( count  ,  1 :  (Acof  +  Bcof+1) )  =  [ Id (n : -1 : (n-Bcof  +  1) )  -l*(Vg. 

((  n-1) : -1: (n-Acof)))  Vg(n)]; 
held(count) .Ms=Ms; 

65  count = count + 1 ; 

end  ; 

"/.solve  the  system  of  equations 
Ns=rref (Ms ) ; 

70 

"/.extract  the  coefficients 
B  =  Ns(l:Bcof  , (Acof+Bcof+1))  ; 

A=[Ns((Bcof+l)  :(Bcof  +  Acof)  ,  (Acof  +  Bcof+1))  ’  ]  ’  ; 

75  "/.regenerate  the  test  input  from  the  linear  difference  ... 

equat ions 

y2= [Vg ( 1 : Bcof )  zeros ( 1 , 400) ] ; 
x2=Id ; 
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7.  iterative  forward  solution 
for  n=(l+max(Acof , Bcof  )  ) : (400) ; 

y2(n)=sum(x2(n:-l:n-Bcof+l) . *B ’ ) -sum (y2 (n-1 : -1 : n-Acof ) . .  . 
.  *  A  ’  )  ; 

end  ; 

st art = length ([Vg(l:Bcof)  zeros(l,floor(max(Acof ,Bcof)+min(... 
Bcof ,Acof)/2))]) ; 

7»error  calculations  &  saving  coefficient  information 

output (q)  . error  =  mean (  (y2 ( start  : 400) -Vg ( start  :  400) ).  “2  ); 

output ( q) . A=A ; 

output ( q) . B=B ; 

output (q) . Acof=Acof ; 

output (q) . Bcof=Bcof ; 

q=q+l ; 

disp (num2str (q) ) ; 

end  ; 


end  ; 

disp ([’ estimates  =  ’  num2str (q-1) ] ) ; 

toe  , 

7«  prepared  all  saved  results  for  error  evaluation 
evals=zeros (1 , length (Alist)*length(Blist)) ; 
f or  n=  1  : 1 :  ( q- 1 )  ; 

evals (n)=output (n) . error ; 

end  ; 

7«  search  out  minima  using  built  in  functions 
Eindex=f ind ( evals ==min ( eval s ) ) ; 
est=l ; 

7«  notify  the  user  if  more  than  one  minima  is  found 
if  length (Eindex >1)  ; 

disp([’found  ’  num2str ( length ( Eindex ) )  ’  mse  mins,  displaying 

first1]); 

end  ; 

7«reload  the  coefficients  for  the  minimum  error  set 
A=output(Eindex(est)) .A; 

B=output(Eindex(est)) .B; 

Acof=output (Eindex (est) ) . Acof ; 

Bcof=output (Eindex (est) ) . Bcof ; 
mse=output (Eindex (est) ) . error ; 

7«  regenerate  the  reconstruction  results  for  display 
y2=[zeros (1 ,420)] ; 
x2= Id ; 
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for  n=(l  +  max(Acof  , Bcof  )  )  :  (400)  ; 

y2 (n) =sum (x2(n:-l:n-Bcof+l)  .  *B  ’  ) -sum (y2(n-l:-l:n-Acof)  ,*A’)  ; 

end  ; 
toe  , 

disp(’N0TE:  The  test  input  has  finite  frequency  content.’); 

disp(’the  resulting  impulse  response  cannot  address  high  frequency 
input  s’); 


"/•plotting  section 
figure ( 1 ) ; 

subplot (2 , 1  ,  1 )  ;  plot (Vg ( 1 : 400) )  ;  axis  tight;  title (’ f iltered  ... 
output ’ ) ; 

subplot (2 , 1 , 2)  ;  plot (y2  ( 1 : 400) )  ;  axis  tight;  title ([’ recon  output 
of  order  A=  ’  num2str ( Acof )  ’  and  B=  ’  num2str ( Bcof ) ] ) ; 

figure (2)  ; 
subplot (2,1,1)  ; 

plot (  ( abs ( Vg ( 1 : 400 ) -y2  ( 1 : 400 ) ) )  ,  ’  r - ’  )  ;  axis  tight;  title(’abs  .. 

error  by  point’); 
subplot (2,1,2)  ; 
plot (Vg  (1  :  400) )  ;  hold  all; 
plot ( y2  (  1  :  400 ))  ;  hold  off; 

title (’ overlaid  voltage  &  reconstructed  voltage’); 

"/•save  data  for  further  use 

save  ABvalsknee2  A  B  Acof  Bcof  Vg  y2  x2 


A.  6  Mat  lab®  Code  for  data  simcomperrord.m 

Listing  A. 6;  simcomperror3.m 

“/•this  process  applies  the  forward  model  coefficents  derived  in 
”/0ABhderiveknee2  to  a  postulated  photonic  input  for  comparison  to 
the 

"/•recorded  signal  &  range  estimation  results  using  Estimation  ... 
Theory 

clear  all;  close  all;  clc; 

"/•  1  .  load  data  for  filter  &  data  set 


"/•this  background  file  was  taken  at  night  making  it  very  close  to  . 
the  dark 

"/•current  values  which  would  normally  be  used. 

[D1,H, stop, mark]=flash_read_seqv2( ’d:\seal\forward  Sim\  ’  .[’plywood 
background  no  subtract2 . seq ’ ] ) ; 


63 


[D2,H, stop, mark]=flash_read_seqv2(’ d:\seal\Forward  Sim\  ’  ,  [’plywood., 
no  subtract  .  seq  ’  ])  ; 

15  '/.initialize  zeros  matrixes  for  target  and  background 
Dtarget  =  zeros (128 , 128 ,20)  ; 

Dbg  =  zeros (128 ,128 ,20)  ; 

°/0 averaging  out  random  noise 
20  for  n=l : 100 ; 

Dbg=Dbg+Dl (n) . data; 

Dtarget=Dtarget+D2 (n) .data; 

end  ; 

25  Dbg=Dbg . /100 ; 

Dtarget  =  Dtarget  .  / 1 0 0 ; 

‘/manually  re-ordering  the  pixels 
Dbg2  =  cat (3 , Dbg ( :  ,  :  ,  3 : 20)  , Dbg ( :  ,  :  ,  1) )  ; 

30  Dtarget2  =  cat(3,Dtarget(:  ,:  ,3:20)  ,Dtarget(:  ,:  ,1)); 

clear  D1  D2  Dbg  Dtarget  H  mark  n  stop; 

7,2.  Find  the  minimum  MSE  for  the  predicted  pulse  shape  estimating  .. 
a  perfect 

35  7, bias  for  first  surface 

%******************************************************************* 


H=210 ; 
hi =105 ; 

40  h2  =  1 04 ; 

firstdata  =  squeeze (Dtarget2(80,64,  :)-Dbg2(80,64,  :))  ’  ; 

MSE  =  zeros (20*180 ,4)  ; 

45  count  =  1  ; 

7,  iteratively  search  through  peak  amplitude  &  pulse  center  ... 
locations  and 

7,calculate  error  between  signal  and  estimate. 

50  for  n= 1 : 20 ; 

for  m=l : 180 ; 

7, bias  estimation  using  all  values  outside  the  pulse  shape 
bias=mean( [firstdata(l: floor ((105+m-hl+10) /30) )  firstdata(. . 
floor ((105  +  m  +  hl  +  10) /30)  :  19) ] )  ; 

55 

7,peak  scaling  solution 
peak=1000+n*10 ; 

7,  scaling  hanning  window  as  pulse  shape 
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pulse=peak*hanning (H)  ; 

60 

'/.building  the  19  point  (currently  upsampled)  estimate 
firstreturn  =  bias*ones (1  ,  19*30)  ; 

f irstreturn (  ( 105+m-hl ) : ( 105+m+h2 )  ) =(pulse ’ ) +bias ; 
f irstdown  =  downsample ( f irstreturn  ,30)  ; 

65 

'/calculate  error  on  downsampled  estimate  compared  to  ... 

recorded 
7.  signal 

MSE ( count  ,  1 ) =n ; 

MSE(count ,2)=m; 

70  MSE ( count  , 3) =mean (  (  f ir st down - f ir st dat a )  .  ~ 2  ); 

MSE(count ,4)=bias; 
count  =  count  + 1 ; 

end  ; 

75  end  ; 

70find  the  minimum  MSE  value 

index  =  find(MSE(:  ,3)==min(MSE(:  ,3))); 

80  "/.pull  out  these  values  for  comparison  and  display  with  the  ... 
original  signal 
n  =  MSE(index(l)  ,1); 
m  =  MSE(index(l)  ,2); 
f irstMSE  =  MSE ( index ( 1 )  ,3)  , 
bias  =  MSE ( index  (  1 )  ,4)  ; 

85 

70  generate  the  estimate 
peak=1000+n*10; 
f ir stpeak  =  peak  , 

90  pulse=peak*hanning (H) ; 

firstreturn  =  bias*ones  (1  ,  19*30)  ; 

f irstreturn (  ( 1 05+m-hl ) : ( 1 05+m+h2 )  )=( pulse ’) +bias ; 

firstcenter=(105+m+30) /30 , 

95 

f irstdown  =  downsample ( f irstreturn  ,30)  ; 

/.plot  the  estimate  and  the  recorded  signal 
figure ( 1 )  ; 

100  subplot (2 , 2 , 1 )  ; 

plot ( f irstdata  ,  1  k :  ’ )  ;  hold  on; 
plot ( f irstdown  ,  1  k- ’  )  ;  hold  off; 
axis  (  [1  19  -200  1500]  )  ; 
title(’first  surface  MSE  +  bias’); 

105  legend ( ’ data ’  ,  ’ sim  ’  )  ; 

7.2 .  Find  the  minimum  MSE  for  the  predicted  pulse  shape  estimating 
a  perfect 
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70bias  for  the  second  surface.  The  code  runs  the  same  except  for  ... 
scaling 


110 

seconddata=squeeze (Dtarget2  (59 ,67  ,  : ) -Dbg2  (59 ,67  ,  : ) )  1  ; 

MSE  =  zeros (20*180 ,4)  ; 
count  =  1  ; 

115 

for  n  =  1  :  2  0  ;  7, 2  0 

for  m=l:  180;  7.20 

bias  =  mean ( [firstdata (1 : floor ( (180  +  m-hl  +  10) /30) )  firstdata(  .  .  . 
floor ((180  +  m  +  hl  +  10) /30)  :  1 9 ) ] )  ; 

120 

peak=1750+n*10  ; 
pulse=peak*hanning (H) ; 
secondreturn=bias*ones(l  ,19*30)  ; 

125  secondreturn  (  ( 180  +  m-hl )  :  ( 180  +  m  +  h2 )  )=( pulse  ’) +bias ; 

seconddown  =  downsample ( secondreturn  ,30)  ; 

MSE ( count , 1 ) =n ; 

MSE(count ,2)=m; 

130  MSE ( count  , 3) =mean (  ( seconddown - seconddata)  .  ~2  ); 

MSE ( count ,4) =bias ; 
count  =  count  + 1 ; 

end  ; 

135  end ; 

index  =  find(MSE(:  ,3)==min(MSE(:  ,3))); 
n  =  MSE(index(l)  ,1); 
m  =  MSE(index(l)  ,2); 

140  secondMSE=MSE ( index ( 1 ) ,3) , 
bias  =  MSE ( index ( 1 )  ,4)  ; 

peak=1750+n*10; 

145  secondpeak  =  peak  , 

pulse=peak*hanning (H) ; 

secondreturn  =  bias*ones (1  ,  19*30)  ; 

secondreturn  (  ( 1 80  +  m-hl )  :  ( 1 80  +  m  +  h2 )  )=( pulse  ’) +bias ; 

150  se condcent er = ( 180+m+30 ) /30 , 

seconddown  =  downsample (secondreturn  ,30)  ; 

figure ( 1 )  ; 
subplot (2,2,2)  ; 

155  plot ( seconddata  ,  1 k :  ’ )  ;  hold  on; 
plot (seconddown  ,  ’k-’)  ;  hold  off; 
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axis  (  [1  19  -200  2000]  )  ; 

title (’ second  surface  MSE  +  bias’); 

legend (’data’  ,  ’  sim’)  ; 

drawnow ; 

7,3 .  now  use  the  pulse  model  system  for  the  gain  point  set  at  the  ... 
knee 

70both  pulse  positions  must  be  considered  to  provide  the  proper  ... 
total 

7«  photonic  demand 


load  ABvalsknee2 ;  clear  y2  Vg  h; 
p=5 . 48 ; 

Vref=XX;  Zvalues  suppressed  for  public  release 
Vbd  =  XX  ;  y.values  suppressed  for  public  release 

bias  =  (4*  1 0  ”  5 )  ;  7.A  value  has  to  be  guessed  for  the  photonic  input 

so  that 

7»gain  variations  outside  the  pulse  can  be  observed 


clear  MSE ; 

MSE  =  zeros (3*3*3*3 ,6)  ; 
count  =  1  ; 
tic  ; 

7«The  values  for  a,b,c,d  can  be  set  closer  to  the  calculated  ... 
results  to 

70display  the  output  quickly  a  =  l:l:3,  b=4:l:6,  c=28:l:30,  d... 
=34: 1:36; 


for  a  =  l : 1 : 20  ; 

for  b  =  1 : 1 : 20; 


for  c  =  1 : 1 : 60  ; 


for  d  =  1 : 1 : 6 0  ; 

70 photonic  pulse  record  building  the  peak  values 
are  based 

70on  quick  radiometric  calculations  ,  but  their  . 
values  are 

70arbitary  as  long  as  they  are  large  enough  to  . 
excite  the 

70gain  variation  process.  The  rescaling  of  the 
algorithm 

% will  match  up  to  the  recorded  results  if  its  . 
relatively  close. 


67 


200 


205 


210 


215 


220 


225 


230 


235 


peakl =160000 ; 

pulse  1  =  peakl * hanning (H )  ’  ; 

peak2  =  160000 ; 
pulse2=peak2*hanning (H) 

photons tempi  =  zeros (1  ,  19*30)  ; 
phot onst emp2  =  zeros (1  ,  19*30)  ; 

phot onst emp 1 (  ( 1 05+ c -hi ) : ( 105+ c+h2 )  )=pulsel; 
phot onst emp2 (  ( 180+d-hl ) : ( 180+d+h2 )  )=pulse2; 

7  compiling  the  photonic  record 

photons  =  [  [zeros (1  ,  120)  downsample ( phot  onst emp 1,30) 
zeros (1 ,41)]  ; 

[zeros  (  1  ,  120)  downsample (photonstemp2  ,  30)  ... 

zeros (1,41)]]; 

"/initialize  the  gain  value  -  only  one  for  the  ... 
array 

M=  [  ( 1 . / (  l-((Vref -0)  ./Vbd)  ,~p  )).* ones  ( 1  ,  180 )]  ; 

"/initialize  a  Vg  &  ID  value  for  each  pixel 
Vg  =  zeros (1  ,180)  ; 

Id  =  zeros (1  ,180)  ; 

"/iteratively  calculate  the  gain  value  due  to  ... 

photonic 
"/interaction  . 

for  n=(l  +  max(Acof  , Bcof ) )  :  1  :  (180)  ; 

Id (n) =sum (M(n)  . *photons ( :  ,n)  ,1)  ; 

Vg  ( n )  =  (  s um  (  I  d  (  n  :  - 1  :  n  -  B  c  o  f  + 1 )  .  *  B  ’  )  -  s um  (  V g  ( n  .  .  . 
-1: -l:n-Acof)  .  *A  ’  ) )  ; 

M(n+l)=(l./(  1 -(( abs (Vref +Vg (n) )). /Vbd) . ~p  )). 
;  °/0M  is  based  on  Vbias 

end  ; 

"/before  adjudicating  the  gain,  we’ll  add  the  ... 

background  level  so  it  gets 
"/modified  assuming  an  otherwise  flat  BG 
for  n= 1 : 2 ; 

Pcount (n, :)=(photons(n, :)+bias) ,*M(1:180) ; 

end  ; 

"/rescale  these  photcounts  to  digital  counts  this  .. 
section 

"/also  helps  reduce  the  required  accuracy  of  the  .  .  . 
photonic 

"/peak  value  estimates 
Mval =(!./(  1 -( (Vref -0) . /Vbd) . ~p  )); 
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f irstdown  =  Pcount  (1  ,  120: 138) -bias* Mval ; 


firstdown=(1160+a*10) *firstdown/(peakl*Mval) ; 

seconddown  =  Pcount (2 , 120 : 138) -bias  +  Mval ; 
seconddown  =  (1850  +  b*10) *  seconddown / (peak2*Mval)  ; 

245 

‘/.calculate  and  store  the  MSEs  for  each  estimated 
pulse 

’/.position  after  applying  the  gain  variation  ... 
process 

MSE(count ,l)=a; 

MSE(count ,2)=b; 

250  MSE ( count , 3) =c ; 

MSE(count ,4)=d; 

MSE(count  ,5)=mean(  ( f ir stdat a - f ir st down )  ,~2  )  ; 
MSE( count  ,6)=mean(  (seconddata-seconddown)  .  ~  2  )  ; 

255  count  =  count  + 1  ; 


end  ; 


260 


end  ; 


end  ; 


end  ; 

265  toe , 

70f  ind  the  minimum  MSEs  for  the  gain  variat  i  onmodel 
indexl=find(  MSE(:,5)==min(MSE(:,5))  ); 
index2=find(  MSE(:,6)==min(MSE(:,6))  ); 

270 

index  1  =  index  1  (1)  ; 
index2  =  index2  (  1 )  ; 

'/.combined  best  values  to  try  to  find  a  common  result 
275  "/.plotting  the  minimums  for  the  second  best  fit 

a  =  MSE(indexl  ,1)  ; 
b  =  MSE ( index2  ,2)  ; 
c  =  MSE(indexl  ,3)  ; 

280  d  =  MSE ( index2  ,4)  ; 

7.  regenerate  the  results  for  display 

peakl =160000 ; 

pulse l=peakl * hanning (H)  ’  ; 

285 


peak2=160000 ; 

pulse2  =  peak2*hanning (H)  ’  ; 

290  phot onst emp 1 =zer os ( 1  ,  1 9*30)  ; 
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photonstemp2=zeros (1  ,19*30)  ; 

phot onst emp 1  (  ( 1 05  + c -hi )  :  ( 1 05+ c  +  h2 )  )=pulsel; 

phot onst emp2 (  ( 1 80+d-hl ) : ( 1 80+d+h2 )  )=pulse2; 

295  f ir st cent er 2  =  ( 105+ c +  60 ) /30  , 
secondcenter2=(180+d+60) /30 , 

photons  =  [ [zeros  (1  ,120)  downsample ( phot onst emp 1  ,30)  zeros  (1,41)]  ; 
[zeros  (  1  ,  120)  downsample ( phot onst emp2  , 30)  zer os ( 1 , 41 ) ] ]  ; 

300 

7,  initialize  the  gain  value  -  only  one  for  the  array 
M=  [  ( 1 . / (  1  - ( ( Vref -0)  . /Vbd)  . ~p  )).* ones ( 1  ,  180) ]  ; 

'/.initialize  a  Vg  &  ID  value  for  each  pixel  140 
305  Vg  =  zeros ( 1  ,  180)  ; 

Id=zeros (1  ,180)  ; 

for  n=(l  +  max(Acof  , Bcof ) )  :  1  :  (180)  ; 

Id (n) =  sum (M(n)  .*photons(:  ,n)  ,1)  ; 

310 

Vg (n)  =  ( sum (Id(n:-l:n-Bcof  +  l)  . *B ’ ) - sum ( Vg (n - 1 : - 1 : n- Acof )  .  *  A  ’  )  ) 
M(n  +  l)=(l./(  1  -((  abs  (  Vref  +Vg  (n)  ))./ Vbd  ).~  p  ));  "/.M  is  based  on 
Vbias 

end  ; 

315  ’/.before  adjudicating  the  gain,  we’ll  add  the  background  level  so 
it  gets 

"/.modified  assuming  an  otherwise  flat  BG 
for  n= 1 : 2 ; 

Pcount (n  ,  :)=(photons(n,  :)+bias)  ,*M(1:180)  ; 

end  ; 

320 


“/.rescale  these  photcounts  to  digital  counts 
Mval=(l./(  l-( (Vref -0) . /Vbd) . "p  )); 

325 

f ir s t do wn  =  Pcount (1  ,  120 : 138) -bias*Mval ; 
f ir st down  =  (1160+a*10) *firstdown/ (peakl*Mval)  ; 

second down  =  Pcount (2  ,  120 : 138) -bias*Mval ; 

330  seconddown  =  (1850  +  b*10) *seconddown/(peak2*Mval)  ; 

f irstMSE2=mean (  ( f ir stdat a - f ir st down )  ,~2  )  ; 
secondMSE2  =  mean (  ( seconddat a - se conddown )  .  "  2  ); 

335  "/.plot  the  results  for  comparison 
figure ( 1 ) ; 
subplot (2,2,3)  ; 

plot ( f irstdata  ,  ’ k :  ’ )  ;  hold  on; 
plot ( f irstdown  ,  ’  k- ’  )  ;  hold  off; 

340  axis (  [1  19  -200  1500]); 
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legend (’data’  ,  ’  sim  ’  )  ; 
title(’first  surface  MSE  +  model’); 

subplot (2,2,4) ; 

345  plot ( seconddata , ’ k : ’ ) ;  hold  on; 
plot ( seconddown  ,  ’  k- ’  )  ;  hold  off; 
axis  (  [1  19  -200  2000]  )  ; 
legend(’data’  ,  ’sim’)  ; 

title(’first  surface  MSE  +  model’); 

350 

7,  save  the  results 
save  simcompresults3 ; 

print  -deps2  simcomperror3.eps; 
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